Transformations.h

Go to the documentation of this file.

#ifndef   math_Transformations_H__
#define   math_Transformations_H__

#include "Transformation.h"
#include "Matrix.h"
#include "utility.h"
#include "../Self.h"

#include <cstdlib>

namespace meow {

template<class Scalar>
class BallProjection: public Transformation<Scalar> {
private:
  struct Myself {
    size_t dimension_;
    Scalar    radius_;

    Myself(size_t d): dimension_(1), radius_(1) {
    }

    Myself(size_t d, Scalar const& r): dimension_(d), radius_(r) {
    }

    Myself(Myself const& m): dimension_(m.dimension_), radius_(m.radius_) {
    }
  };

  Self<Myself> const self;
public:
  BallProjection(BallProjection const& b): Transformation<Scalar>(b),
  self(b.self, Self<Myself>::COPY_FROM) {
  }

  BallProjection(size_t d): Transformation<Scalar>(d, 1, d, 1, 1),
  self(Myself(d)) {
    radius(1);
  }

  BallProjection(size_t d, Scalar const& r): Transformation<Scalar>(d,1,d,1,1),
  self(Myself(d, r)) {
    radius(r);
  }

  BallProjection& copyFrom(BallProjection const& b) {
    Transformation<Scalar>::copyFrom(b);
    copyFrom(b);
    return *this;
  }

  BallProjection& referenceFrom(BallProjection const& b) {
    Transformation<Scalar>::referenceFrom(b);
    referenceFrom(b);
    return *this;
  }

  Scalar parameter(size_t i) const {
    return radius();
  }

  Scalar parameter(size_t i, Scalar const& s) {
    return radius(s);
  }

  Scalar radius() const {
    return self->radius_;
  }

  Scalar radius(Scalar const& r) {
    self()->radius_ = r;
    return radius();
  }

  size_t dimension() const {
    return self->dimension_;
  }


  Matrix<Scalar> transformate(Matrix<Scalar> const& x) const {
    Matrix<Scalar> ret(x);
    for (size_t c = 0, C = ret.cols(); c < C; c++) {
      Scalar sum(0);
      for (size_t i = 0; i < self->dimension_; i++) {
        sum = sum + squ(ret(i, c));
      }
      Scalar len(sqrt(double(sum)));
      for (size_t i = 0; i < self->dimension_; i++) {
        ret(i, c, ret(i, c) * radius() / len);
      }
    }
    return ret;
  }

  Matrix<Scalar> jacobian(Matrix<Scalar> const& x) const {
    Scalar sum(0);
    for(size_t i = 0, I = dimension(); i < I; ++i)
      sum = sum + squ(x(i, 0));
    Scalar len(sqrt(double(sum)));
    Matrix<Scalar> ret(dimension(), dimension(), Scalar(0.0));
    for(size_t i = 0, I = dimension(); i < I; ++i)
      for(size_t j = 0; j < I; ++j)
        if (i == j) {
          ret(i, j, radius() * (squ(len) - squ(x(i, 0))) / cub(len));
        }
        else {
          ret(i, j, radius() * (-x(i, 0) * x(j, 0) / cub(len)));
        }
    return ret;
  }

  Matrix<Scalar> jacobian(Matrix<Scalar> const& x, size_t i) const {
    Matrix<Scalar> ret(dimension(), 1, Scalar(0.0));
    Scalar sum(0);
    for(size_t i = 0, I = dimension(); i < I; i++) {
      sum = sum + squ(x(i, 0));
    }
    return ret / Scalar(sqrt(double(sum)));
  }

  BallProjection& operator=(BallProjection const& b) {
    return copyFrom(b);
  }

  Matrix<Scalar> operator()(Matrix<Scalar> const& v) const {
    return transformate(v);
  }
};


template<class Scalar>
class PhotoProjection: public Transformation<Scalar> {
private:
  struct Myself {
    Scalar     focal_;
    size_t dimension_;

    Myself() {
    }

    Myself(size_t d, Scalar f): focal_(f), dimension_(d) {
    }

    Myself(Myself const& b): focal_(b.focal_), dimension_(b.dimension_) {
    }

    ~Myself() {
    }
  };

  Self<Myself> const self;
public:
  PhotoProjection(size_t dimension):
  Transformation<Scalar>(dimension, 1, dimension, 1, 1),
  self(Myself(dimension, 1)) {
  }

  PhotoProjection(size_t dimension, Scalar const& f):
  Transformation<Scalar>(dimension, 1, dimension, 1, 1),
  self(Myself(dimension, f)) {
  }

  PhotoProjection(PhotoProjection const& p): Transformation<Scalar>(p),
  self(p.self, Self<Myself>::COPY_FROM) {
  }

  PhotoProjection& copyFrom(PhotoProjection const& b) {
    Transformation<Scalar>::copyFrom(b);
    self().copyFrom(b.self);
    return *this;
  }

  PhotoProjection& referenceFrom(PhotoProjection const& b) {
    Transformation<Scalar>::referenceFrom(b);
    self().referenceFrom(b.self);
    return *this;
  }

  Scalar parameter(size_t i) const {
    return focal();
  }

  Scalar parameter(size_t i, Scalar const& s){
    return focal(s);
  }

  Scalar focal() const {
    return self->focal_;
  }

  Scalar focal(Scalar const& f){
    self()->focal_ = f;
    return focal();
  }

  size_t dimension() const {
    return self->dimension_;
  }

  Matrix<Scalar> transformate(Matrix<Scalar> const& x) const {
    Matrix<Scalar> ret(x);
    for (size_t c = 0, C = ret.cols(); c < C; c++) {
      for (size_t i = 0, I = dimension(); i < I; ++i) {
        ret(i, c, -ret(i, c) * focal() / ret(I - 1, c));
      }
    }
    return ret;
  }

  Matrix<Scalar> jacobian(Matrix<Scalar> const& x) const{
    Matrix<Scalar> ret(dimension(), dimension(), Scalar(0.0));
    for(ssize_t i = 0, I = (ssize_t)dimension() - 1; i < I; i++){
      ret(i,               i, -focal() /     x(I, 0) );
      ret(i, dimension() - 1,  focal() / squ(x(I, 0)));
    }
    return ret;
  }

  Matrix<Scalar> jacobian(Matrix<Scalar> const& x, size_t i) const{
    Matrix<Scalar> ret(dimension(), 1, Scalar(0.0));
    for(size_t i = 0, I = dimension(); i < I; ++i) {
      ret(i, 0, -x(i, 0) / x(I - 1, 0));
    }
    return ret;
  }

  PhotoProjection& operator=(PhotoProjection const& b) {
    return copyFrom(b);
  }

  Matrix<Scalar> operator()(Matrix<Scalar> const& v) const {
    return transformate(v);
  }
};

} // meow

#endif // Transformations_H__