TextureTransform.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 <Base/Common/BIASpragmaStart.hh>
#include <Image/TextureTransform.hh>

int BIAS::TextureTransform::
TextureJacobianForward(const HomgPoint2D& src, 
                       Matrix2x2<double>& Jac) 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 sink_left, sink_right, sink_up, sink_down;
  HomgPoint2D loc(src);
  loc[0] -= 0.5;
  if (MapForward(loc, sink_left)) return -1;
  loc[0] += 1.0;
  if (MapForward(loc, sink_right)) return -1;
  loc[0] -= 0.5;
  loc[1] -= 0.5;
  if (MapForward(loc, sink_up)) return -1;
  loc[1] += 1.0;
  if (MapForward(loc, sink_down)) return -1;

  // difference vectors are rows of jacobian
  Jac[0][0] = sink_right[0] - sink_left[0];
  Jac[1][0] = sink_right[1] - sink_left[1];
  Jac[0][1] = sink_down[0]  - sink_up[0];
  Jac[1][1] = sink_down[1]  - sink_up[1];
  return 0;
}

int BIAS::TextureTransform::
TextureJacobianBackward(const HomgPoint2D& sink, 
                        Matrix2x2<double>& Jac) 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 (MapBackward(locsink, source_left)) return -1;
  locsink[0] += 1.0;
  if (MapBackward(locsink, source_right)) return -1;
  locsink[0] -= 0.5;
  locsink[1] -= 0.5;
  if (MapBackward(locsink, source_up)) return -1;
  locsink[1] += 1.0;
  if (MapBackward(locsink, source_down)) return -1;

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

int BIAS::TextureTransform::
ParameterJacobianForward(Matrix<double>& Jac, const HomgPoint2D& src) {
  // remember current parameters
  Vector<double> p;
  GetParameters(p);
  const unsigned int numvars = p.Size();
  Jac.newsize(2,numvars);
  HomgPoint2D p1, ref;
  // map current point with current parameters
  MapForward(src, ref);

  // vary each entry in p by 2%
  for (unsigned int j=0; j<numvars; j++) {
    Vector<double> pchanged(p);
    double myeps = 0.02 * fabs(p[j]);
    if (myeps<1e-9) myeps = 1e-9;
    pchanged[j] += myeps;
    // set the "new" parameters 
    SetParameters(pchanged);
    // and map
    MapForward(src, p1);
    // compute difference and approximate jacobian entry
    p1 -= ref;
    Jac[0][j] = p1[0]/myeps;
    Jac[1][j] = p1[0]/myeps;
  }
  // reset original parameters
  SetParameters(p);
  return 0;
}

int BIAS::TextureTransform::
ParameterJacobianBackward(Matrix<double>& Jac, const HomgPoint2D& sink) {
  // remember current parameters
  Vector<double> p;
  GetParameters(p);
  const unsigned int numvars = p.Size();
  Jac.newsize(2, numvars);
  HomgPoint2D p1, ref;
  // map current point with current parameters
  MapBackward(sink, ref);

  // vary each entry in p by 2%
  for (unsigned int j=0; j<numvars; j++) {
    Vector<double> pchanged(p);
    double myeps = 0.02 * fabs(p[j]);
    if (myeps<1e-9) myeps = 1e-9;
    pchanged[j] += myeps;
    // set the "new" parameters 
    SetParameters(pchanged);
    // and map
    MapBackward(sink, p1);
    // compute difference and approximate jacobian entry
    p1 -= ref;
    Jac[0][j] = p1[0]/myeps;
    Jac[1][j] = p1[0]/myeps;
  }
  // reset original parameters
  SetParameters(p);
  return 0;
}