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 {
    Scalar    radius_;
    size_t dimension_;

    Myself() {
    }
    ~Myself() {
    }
    Myself& copyFrom(Myself const& b) {
      radius_    = b.radius_;
      dimension_ = b.dimension_;
      return *this;
    }
  };
  
  Self<Myself> const self;
public:
  BallProjection(BallProjection const& b):
  Transformation<Scalar>(b), self(false) {
    copyFrom(b);
  }
  
  BallProjection(size_t d): self(true),
  Transformation<Scalar>(d, 1, d, 1, 1) {
    self()->dimension_ = d;
    radius(1);
  }
  
  BallProjection(size_t d, Scalar const& r):
  Transformation<Scalar>(d, 1, d, 1, 1), self(true) {
    self()->dimension_ = d;
    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() {
    }
    Myself& copyFrom(Myself const& b) {
      focal_     = b.focal_;
      dimension_ = b.dimension_;
      return *this;
    }
  };
  
  Self<Myself> const& self;
public:
  PhotoProjection(size_t dimension): 
  Transformation<Scalar>(dimension, 1, dimension, 1, 1), self(true) {
    self()->dimension_ = dimension;
    focal(1);
  }
  
  PhotoProjection(size_t dimension, Scalar const& f):
  Transformation<Scalar>(dimension, 1, dimension, 1, 1), self(true) {
    self()->dimension_ = dimension;
    focal(f);
  }
  
  PhotoProjection(PhotoProjection const& p):
  Transformation<Scalar>(p), self(false) {
    self().copyFrom(p.self);
  }
  
  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);
  }
};

}

#endif // Transformations_H__