diff options
Diffstat (limited to 'meowpp/math')
| -rw-r--r-- | meowpp/math/!readme.asciidoc | 73 | ||||
| -rw-r--r-- | meowpp/math/LinearTransformation.h | 110 | ||||
| -rw-r--r-- | meowpp/math/LinearTransformations.h | 404 | ||||
| -rw-r--r-- | meowpp/math/Matrix.h | 536 | ||||
| -rw-r--r-- | meowpp/math/Transformation.h | 237 | ||||
| -rw-r--r-- | meowpp/math/Transformations.h | 550 | ||||
| -rw-r--r-- | meowpp/math/Vector.h | 267 | ||||
| -rw-r--r-- | meowpp/math/methods.h | 230 | ||||
| -rw-r--r-- | meowpp/math/utility.h | 157 |
9 files changed, 0 insertions, 2564 deletions
diff --git a/meowpp/math/!readme.asciidoc b/meowpp/math/!readme.asciidoc deleted file mode 100644 index 062b45d..0000000 --- a/meowpp/math/!readme.asciidoc +++ /dev/null @@ -1,73 +0,0 @@ - - -===== utility.h - -數學相關的小 function 雜七雜八的不知道歸類何處 - -.Functions -* noEPS() -* normalize() -* denormalize() -* ratioMapping() -* inRange() -* squ() -* cub() -* average() -* average() -* tAbs() - -.Constants -* PI - -===== Matrix.h - -.Classes -* `meow::Matrix<Entry>` - -===== Vector.h - -實作上將 *Matrix* 重新包裝 - -.Classes -* `meow::Vector<Scalar>` - -===== Transformation.h - -各種轉換的 Base Class, 這裡所謂的 *Transformation* 形式上不一定要是 Linear, -但原則上都是 *input a vector, output a vector* 其中input/output的dimension可以 -不同. - -.Classes -* `meow::Transformation<Scalar>` - -===== Transformations.h - -包含各種 *Non-Linear* transformation - -.Classes -* `meow::BallProjection<Scalar>` -* `meow::PhotoProjection<Scalar>` - -===== LinearTransformation.h - -各種 LinearTransformation 的Base Class, 繼承自 `meow::Transformation` - -.Classes -* `meow::LinearTransformation<Scalar>` - -===== LinearTransformations.h - -各種 *Linear* Transformation - -.Classes -* `meow::Rotation3D<Scalar>` - -===== methods.h - -一些數學方法 - -.Functions -* ransac() -* levenbergMarquardt() - - diff --git a/meowpp/math/LinearTransformation.h b/meowpp/math/LinearTransformation.h deleted file mode 100644 index 0dc3735..0000000 --- a/meowpp/math/LinearTransformation.h +++ /dev/null @@ -1,110 +0,0 @@ -#ifndef math_LinearTransformation_H__ -#define math_LinearTransformation_H__ - -#include "Transformation.h" -#include "Matrix.h" - -#include <cstdlib> - -namespace meow { - -/*! - * @brief A base class for implementing kinds of linear transformations. - * - * Because all linear transformations belong to transformations, - * this class inherit to Transformation. - * - * @author cat_leopard - */ -template<class Scalar> -class LinearTransformation: public Transformation<Scalar> { -private: - Matrix<Scalar> matrix_; -protected: - /*! - * Constructor with input/output size gived - */ - LinearTransformation(size_t inputRows, size_t outputRows, size_t psize): - Transformation<Scalar>(inputRows, 1u, outputRows, 1u, psize), - matrix_(outputRows, inputRows, Scalar(0.0)) { - } - - /*! - * Constructor with input/output size gived and a inital matrix - */ - LinearTransformation(size_t inputRows, size_t outputRows, size_t psize, - Matrix<Scalar> const& m): - Transformation<Scalar>(inputRows, 1u, outputRows, 1u, psize), - matrix_(m) { - } - - /*! - * Constructor with another LinearTransformation - * - * @param [in] b another LinearTransformation - */ - LinearTransformation(LinearTransformation const& b): - Transformation<Scalar>(b), - matrix_(b.matrix_) { - } - - /*! - * @brief Copy settings, matrix from another LinearTransformation - * - * @param [in] b another LinearTransformation - */ - LinearTransformation& copyFrom(LinearTransformation const& b) { - Transformation<Scalar>::copyFrom(b); - matrix_.copyFrom(b.matrix_); - return *this; - } - - /*! - * @brief Reference settings, matrix from another LinearTransformation - * - * @param [in] b another LinearTransformation - */ - LinearTransformation& referenceFrom(LinearTransformation const& b) { - Transformation<Scalar>::referenceFrom(b); - matrix_.referenceFrom(b.matrix_); - return *this; - } - - /*! - * @brief setup the matrix - */ - virtual Matrix<Scalar> const& matrix(Matrix<Scalar> const& m) { - matrix_.copyFrom(m); - return matrix(); - } - -public: - /*! - * Destructor - */ - virtual ~LinearTransformation() { - } - - /*! - * @brief Return the matrix form of this transformation - * - * @return A matrix - */ - virtual Matrix<Scalar> const& matrix() const { - return matrix_; - } - - /*! - * @brief Return the inverse of the matrix form of this transformate - * - * @return A matrix (may be invalid) - */ - virtual Matrix<Scalar> matrixInv() const { - return matrix_.inverse(); - } -}; - - -} // meow - -#endif // math_LinearTransformation_H__ diff --git a/meowpp/math/LinearTransformations.h b/meowpp/math/LinearTransformations.h deleted file mode 100644 index 4bf9a36..0000000 --- a/meowpp/math/LinearTransformations.h +++ /dev/null @@ -1,404 +0,0 @@ -#ifndef math_LinearTransformations_H__ -#define math_LinearTransformations_H__ - -#include "LinearTransformation.h" -#include "Matrix.h" -#include "utility.h" -#include "../Self.h" -#include "../geo/Vectors.h" - -#include <cstdlib> - -namespace meow { - -/*! - * @brief Rotation a point/vector alone an axis with given angle in 3D world. - * - * @author cat_leopard - */ -template<class Scalar> -class Rotation3D: public LinearTransformation<Scalar> { -private: - struct Myself { - Vector3D<Scalar> theta_; - bool need_; - - Myself(): theta_(0, 0, 0), need_(true) { - } - - Myself(Myself const& b): theta_(b.theta_), need_(b.need_) { - } - - ~Myself() { - } - }; - - Self<Myself> const self; - - void calcMatrix() const { - if (self->need_) { - Matrix<Scalar> tmp(3, 3, 0.0); - if (noEPS(self->theta_.length2()) == 0) { - tmp.identitied(); - } - else { - Vector3D<double> axis (self->theta_.normalize()); - double angle(self->theta_.length()); - double cs(cos(angle / 2.0)); - double sn(sin(angle / 2.0)); - - tmp.entry(0, 0, 2*(squ(axis.x())-1.0)*squ(sn) + 1); - tmp.entry(1, 1, 2*(squ(axis.y())-1.0)*squ(sn) + 1); - tmp.entry(2, 2, 2*(squ(axis.z())-1.0)*squ(sn) + 1); - tmp.entry(0, 1, 2*axis.x()*axis.y()*squ(sn) - 2*axis.z()*cs*sn); - tmp.entry(1, 0, 2*axis.y()*axis.x()*squ(sn) + 2*axis.z()*cs*sn); - tmp.entry(0, 2, 2*axis.x()*axis.z()*squ(sn) + 2*axis.y()*cs*sn); - tmp.entry(2, 0, 2*axis.z()*axis.x()*squ(sn) - 2*axis.y()*cs*sn); - tmp.entry(1, 2, 2*axis.y()*axis.z()*squ(sn) - 2*axis.x()*cs*sn); - tmp.entry(2, 1, 2*axis.z()*axis.y()*squ(sn) + 2*axis.x()*cs*sn); - } - ((Rotation3D*)this)->LinearTransformation<Scalar>::matrix(tmp); - self()->need_ = false; - } - } - -public: - /*! - * Constructor with no rotation - */ - Rotation3D(): LinearTransformation<Scalar>(3u, 3u, 3u), self() { - } - - /*! - * Constructor and copy data - */ - Rotation3D(Rotation3D const& b): LinearTransformation<Scalar>(b), - self(b.self, Self<Myself>::COPY_FROM) { - } - - /*! - * Destructor - */ - ~Rotation3D() { - } - - /*! - * @brief Copy data - * - * @param [in] b another Rotation3D class. - * @return \c *this - */ - Rotation3D& copyFrom(Rotation3D const& b) { - LinearTransformation<Scalar>::copyFrom(b); - self().copyFrom(b.self); - return *this; - } - - /*! - * @brief Reference data - * - * @param [in] b another Rotation3D class. - * @return \c *this - */ - Rotation3D& referenceFrom(Rotation3D const& b) { - LinearTransformation<Scalar>::referenceFrom(b); - self().referenceFrom(b.self); - return *this; - } - - /*! - * @brief same as \c theta(i) - */ - Scalar parameter(size_t i) const { - return theta(i); - } - - /*! - * @brief same as \c theta(i, s) - */ - Scalar parameter(size_t i, Scalar const& s) { - return theta(i, s); - } - - /*! - * @brief Get the \c i -th theta - * - * \c i can only be 1, 2 or 3 - * - * @param [in] i index - * @return \c i -th theta - */ - Scalar const& theta(size_t i) const { - return self->theta_(i); - } - - /*! - * @brief Set the \c i -th theta - * - * \c i can only be 1, 2 or 3 - * - * @param [in] i index - * @param [in] s new theta value - * @return \c i -th theta - */ - Scalar const& theta(size_t i, Scalar const& s) { - if (theta(i) != s) { - if (i == 0) self()->theta_.x(s); - else if (i == 1) self()->theta_.y(s); - else if (i == 2) self()->theta_.z(s); - self()->need_ = true; - } - return theta(i); - } - - /*! - * @brief Setting - * - * @param [in] axis axis - * @param [in] angle angle - */ - void axisAngle(Vector<Scalar> const& axis, Scalar const& angle) { - Vector<Scalar> n(axis.normalize()); - for (size_t i = 0; i < 3; i++) { - theta(i, n(i) * angle); - } - } - - /*! - * @brief Concat another rotation transformation - * @param [in] r another rotation transformation - */ - Rotation3D& add(Rotation3D const& r) { - for (size_t i = 0; i < 3; i++) { - theta(i, r.theta(i)); - } - return *this; - } - - /*! - * @brief Do the transformate - - * Assume: - * - The input vector is \f$ (x ,y ,z ) \f$ - * - The output vector is \f$ (x',y',z') \f$ - * - The parameters theta is\f$ \vec{\theta}=(\theta_x,\theta_y,\theta_z) \f$ - * . - * Then we have: - * \f[ - * \left[ \begin{array}{c} x' \\ y' \\ z' \\ \end{array} \right] - * = - * \left[ \begin{array}{ccc} - * 2(n_x^2 - 1) \sin^2\phi + 1 & - * 2n_x n_y \sin^2\phi - 2n_z\cos \phi\sin \phi & - * 2n_x n_z \sin^2\phi + 2n_y\cos \phi\sin \phi \\ - * 2n_y n_x \sin^2\phi + 2n_z\cos \phi\sin \phi & - * 2(n_y^2 - 1) \sin^2\phi + 1 & - * 2n_y n_z \sin^2\phi - 2n_x\cos \phi\sin \phi \\ - * 2n_z n_x \sin^2\phi - 2n_y\cos \phi\sin \phi & - * 2n_z n_y \sin^2\phi + 2n_x\cos \phi\sin \phi & - * 2(n_z^2 - 1) \sin^2\phi + 1 \\ - * \end{array} \right] - * \left[ \begin{array}{c} x \\ y \\ z \\ \end{array} \right] - * \f] - * Where: - * - \f$ \phi \f$ is the helf of length of \f$ \vec{\theta} \f$ , - * which means \f$ \phi = \frac{\left|\vec{\theta}\right|}{2} - * = \frac{1}{2}\sqrt{\theta_x^2 + \theta_y^2 + \theta_z^2} \f$ - * - \f$ \vec{n} \f$ is the normalized form of \f$ \vec{\theta} \f$ , - * which means \f$ \vec{n} = (n_x,n_y,n_z) = \vec{\theta} / 2\phi \f$ - * - * @param [in] x the input vector - * @return the output matrix - */ - Matrix<Scalar> transformate(Matrix<Scalar> const& x) const { - calcMatrix(); - return LinearTransformation<Scalar>::matrix() * x; - } - - /*! - * @brief Return the jacobian matrix (derivate by the input vector) - * of this transformate - * - * The matrix we return is: - * \f[ - * \left[ \begin{array}{ccc} - * 2(n_x^2 - 1) \sin^2\phi + 1 & - * 2n_x n_y \sin^2\phi - 2n_z\cos \phi\sin \phi & - * 2n_x n_z \sin^2\phi + 2n_y\cos \phi\sin \phi \\ - * 2n_y n_x \sin^2\phi + 2n_z\cos \phi\sin \phi & - * 2(n_y^2 - 1) \sin^2\phi + 1 & - * 2n_y n_z \sin^2\phi - 2n_x\cos \phi\sin \phi \\ - * 2n_z n_x \sin^2\phi - 2n_y\cos \phi\sin \phi & - * 2n_z n_y \sin^2\phi + 2n_x\cos \phi\sin \phi & - * 2(n_z^2 - 1) \sin^2\phi + 1 \\ - * \end{array} \right] - * \f] - * Where the definition of \f$ \vec{n} \f$ and \f$ \phi \f$ - * is the same as the definition in the description of - * the method \b transformate() . - * - * @param [in] x the input vector (in this case it is a useless parameter) - * @return a matrix - */ - Matrix<Scalar> jacobian(Matrix<Scalar> const& x) const { - calcMatrix(); - return LinearTransformation<Scalar>::matrix(); - } - - /*! - * @brief Return the jacobian matrix of this transformate - * - * Here we need to discussion in three case: - * - \a i = 0, derivate by the x axis of the vector theta - * \f[ - * \left[ \begin{array}{ccc} - * 0 & 0 & 0 \\ - * 0 & 0 & -1 \\ - * 0 & 1 & 0 \\ - * \end{array} \right] - * \left[ \begin{array}{ccc} - * 2(n_x^2 - 1) \sin^2\phi + 1 & - * 2n_x n_y \sin^2\phi - 2n_z\cos \phi\sin \phi & - * 2n_x n_z \sin^2\phi + 2n_y\cos \phi\sin \phi \\ - * 2n_y n_x \sin^2\phi + 2n_z\cos \phi\sin \phi & - * 2(n_y^2 - 1) \sin^2\phi + 1 & - * 2n_y n_z \sin^2\phi - 2n_x\cos \phi\sin \phi \\ - * 2n_z n_x \sin^2\phi - 2n_y\cos \phi\sin \phi & - * 2n_z n_y \sin^2\phi + 2n_x\cos \phi\sin \phi & - * 2(n_z^2 - 1) \sin^2\phi + 1 \\ - * \end{array} \right] - * \left[ \begin{array}{c} x \\ y \\ z \\ \end{array} \right] - * \f] - * - \a i = 1, derivate by the y axis of the vector theta - * \f[ - * \left[ \begin{array}{ccc} - * 0 & 0 & 1 \\ - * 0 & 0 & 0 \\ - * -1 & 0 & 0 \\ - * \end{array} \right] - * \left[ \begin{array}{ccc} - * 2(n_x^2 - 1) \sin^2\phi + 1 & - * 2n_x n_y \sin^2\phi - 2n_z\cos \phi\sin \phi & - * 2n_x n_z \sin^2\phi + 2n_y\cos \phi\sin \phi \\ - * 2n_y n_x \sin^2\phi + 2n_z\cos \phi\sin \phi & - * 2(n_y^2 - 1) \sin^2\phi + 1 & - * 2n_y n_z \sin^2\phi - 2n_x\cos \phi\sin \phi \\ - * 2n_z n_x \sin^2\phi - 2n_y\cos \phi\sin \phi & - * 2n_z n_y \sin^2\phi + 2n_x\cos \phi\sin \phi & - * 2(n_z^2 - 1) \sin^2\phi + 1 \\ - * \end{array} \right] - * \left[ \begin{array}{c} x \\ y \\ z \\ \end{array} \right] - * \f] - * - \a i = 2, derivate by the z axis of the vector theta - * \f[ - * \left[ \begin{array}{ccc} - * 0 & -1 & 0 \\ - * 1 & 0 & 0 \\ - * 0 & 0 & 0 \\ - * \end{array} \right] - * \left[ \begin{array}{ccc} - * 2(n_x^2 - 1) \sin^2\phi + 1 & - * 2n_x n_y \sin^2\phi - 2n_z\cos \phi\sin \phi & - * 2n_x n_z \sin^2\phi + 2n_y\cos \phi\sin \phi \\ - * 2n_y n_x \sin^2\phi + 2n_z\cos \phi\sin \phi & - * 2(n_y^2 - 1) \sin^2\phi + 1 & - * 2n_y n_z \sin^2\phi - 2n_x\cos \phi\sin \phi \\ - * 2n_z n_x \sin^2\phi - 2n_y\cos \phi\sin \phi & - * 2n_z n_y \sin^2\phi + 2n_x\cos \phi\sin \phi & - * 2(n_z^2 - 1) \sin^2\phi + 1 \\ - * \end{array} \right] - * \left[ \begin{array}{c} x \\ y \\ z \\ \end{array} \right] - * \f] - * . - * Where \f$ (x,y,z) \f$ is the input vector, \f$ \vec{n}, \phi \f$ is the - * same one in the description of \b transformate(). - * - * @param [in] x the input vector - * @param [in] i the index of the parameters(theta) to dervite - * @return a matrix - */ - Matrix<Scalar> jacobian(Matrix<Scalar> const& x, size_t i) const { - calcMatrix(); - Matrix<Scalar> mid(3u, 3u, Scalar(0.0)); - if (i == 0) { - mid.entry(1, 2, Scalar(-1.0)); - mid.entry(2, 1, Scalar( 1.0)); - } - else if(i == 1) { - mid.entry(0, 2, Scalar( 1.0)); - mid.entry(2, 0, Scalar(-1.0)); - } - else { - mid.entry(0, 1, Scalar(-1.0)); - mid.entry(1, 0, Scalar( 1.0)); - } - return mid * LinearTransformation<Scalar>::matrix() * x; - } - - /*! - * @brief Do the inverse transformate - * - * @param [in] x the input vector - * @return the output vector - */ - Matrix<Scalar> transformateInv(Matrix<Scalar> const& x) const { - return matrixInv() * x; - } - - /*! - * @brief Return the jacobian matrix of the inverse form of this transformate - * - * @param [in] x the input vector - * @return a matrix - */ - Matrix<Scalar> jacobianInv(Matrix<Scalar> const& x) const { - return matrixInv(); - } - - /*! - * @brief Return the jacobian matrix of the inverse form of this transformate - * - * @param [in] x the input vector - * @param [in] i the index of the parameters(theta) to dervite - * @return a matrix - */ - Matrix<Scalar> jacobianInv(Matrix<Scalar> const& x, size_t i) const { - calcMatrix(); - Matrix<Scalar> mid(3u, 3u, Scalar(0.0)); - if (i == 0) { - mid.entry(1, 2, Scalar(-1.0)); - mid.entry(2, 1, Scalar( 1.0)); - } - else if(i == 1) { - mid.entry(0, 2, Scalar( 1.0)); - mid.entry(2, 0, Scalar(-1.0)); - } - else { - mid.entry(0, 1, Scalar(-1.0)); - mid.entry(1, 0, Scalar( 1.0)); - } - return matrixInv() * mid.transpose() * x; - return (-mid) * matrixInv() * x; - } - - /*! - * @brief Return the inverse matrix - * - * In this case, the inverse matrix is equal to the transpose of the matrix - * - * @return a matrix - */ - Matrix<Scalar> matrixInv() const { - calcMatrix(); - return LinearTransformation<Scalar>::matrix().transpose(); - } - - //! @brief same as \c copyFrom(b) - Rotation3D& operator=(Rotation3D const& b) { - return copyFrom(b); - } -}; - -} // meow - -#endif // math_LinearTransformations_H__ diff --git a/meowpp/math/Matrix.h b/meowpp/math/Matrix.h deleted file mode 100644 index a2f45aa..0000000 --- a/meowpp/math/Matrix.h +++ /dev/null @@ -1,536 +0,0 @@ -#ifndef math_Matrix_H__ -#define math_Matrix_H__ - -#include "../Self.h" -#include "../math/utility.h" - -#include <vector> -#include <algorithm> - -#include <cstdlib> - -namespace meow { -/*! - * @brief \b matrix - * - * @author cat_leopard - */ -template<class Entry> -class Matrix { -public: - typedef typename std::vector<Entry>::reference EntryRef ; - typedef typename std::vector<Entry>::const_reference EntryRefK; -private: - struct Myself { - size_t rows_; - size_t cols_; - std::vector<Entry> entries_; - - Myself(): - rows_(0), cols_(0), entries_(0) { - } - - Myself(Myself const& b): - rows_(b.rows_), cols_(b.cols_), entries_(b.entries_) { - } - - Myself(size_t r, size_t c, Entry const& e): - rows_(r), cols_(c), entries_(r * c, e) { - } - - ~Myself() { - } - - size_t index(size_t r, size_t c) const { - return r * cols_ + c; - } - - void realSize() { - std::vector<Entry> tmp(entries_); - entries_.swap(tmp); - } - }; - - Self<Myself> const self; -public: - /*! - * @brief constructor - * - * Create an empty matrix with size \b 0x0. - * In other world, create an \b invalid matrix - */ - Matrix(): self() { } - - /*! - * @brief constructor - * - * Copy data from another one - * - * @param [in] m another matrix - */ - Matrix(Matrix const& m): self(m.self, Self<Myself>::COPY_FROM) { - } - - /*! - * @brief constructor - * - * Create an \a r x \a c matrix with all entry be \a e - * - * @param [in] r number of rows - * @param [in] c number of columns - * @param [in] e inital entry - */ - Matrix(size_t r, size_t c, Entry const& e): self(Myself(r, c, e)) { - } - - //! @brief destructor - ~Matrix() { } - - /*! - * @brief copy - * - * Copy data from another matrix - * - * @param [in] m matrix - * @return *this - */ - Matrix& copyFrom(Matrix const& m) { - self().copyFrom(m.self); - return *this; - } - - /*! - * @brief reference - * - * Reference itself to another matrix - * - * @param [in] m matrix - * @return *this - */ - Matrix& referenceFrom(Matrix const& m) { - self().referenceFrom(m.self); - return *this; - } - - //! @brief reset the size of the matrix to \a r x \a c with entry all be \a e - void reset(size_t r, size_t c, Entry const& e) { - self()->rows_ = r; - self()->cols_ = c; - self()->entries_.clear(); - self()->entries_.resize(r * c, e); - } - - //! @brief Return whether it is a \b valid matrix - bool valid() const { - return (rows() > 0 && cols() > 0); - } - - //! @brief Return number of rows - size_t rows() const { - return self->rows_; - } - - //! @brief Return number of cols - size_t cols() const { - return self->cols_; - } - - //! @brief Return number of rows times number of cols - size_t size() const { - return rows() * cols(); - } - - /*! - * @brief resize the matrix such that number of rows become \a r. - * - * New created entry will be \a e - * - * @param [in] r new number of rows - * @param [in] e inital entry - * @return new number of rows - */ - size_t rows(size_t r, Entry const& e) { - if (r != rows()) { - self()->entries_.resize(r * cols(), e); - self()->rows_ = r; - } - return rows(); - } - - /*! - * @brief resize the matrix such that number of cols become \a c - * - * New created entry will be \a e - * - * @param [in] c new number of columns - * @param [in] e inital entry - * @return new number of columns - */ - size_t cols(size_t c, Entry const& e) { - if (c != cols()) { - Self<Myself> const old(self, Self<Myself>::COPY_FROM); - self()->entries_.resize(rows() * c); - self()->cols_ = c; - for (size_t i = 0, I = rows(); i < I; i++) { - size_t j, J1 = std::min(old->cols_, cols()), J2 = cols(); - for (j = 0; j < J1; j++) - self()->entries_[self->index(i, j)] = old->entries_[old->index(i, j)]; - for (j = J1; j < J2; j++) - self()->entries_[self->index(i, j)] = e; - } - } - return cols(); - } - - /*! - * @brief resize - * - * Resize to \a r x \a c, with new created entry be \a e - * - * @param [in] r number of rows - * @param [in] c number of rows - * @param [in] e inital entry - * @return \a r * \a c - */ - size_t size(size_t r, size_t c, Entry const& e) { - cols(c, e); - rows(r, e); - return rows() * cols(); - } - - /*! - * @brief free the memory - */ - void clear() { - self()->rows_ = 0; - self()->cols_ = 0; - self()->entries_.clear(); - self()->realSize(); - } - - //! @brief Access the entry at \a r x \a c - Entry entry(size_t r, size_t c) const { - return self->entries_[self->index(r, c)]; - } - - //! @brief Change the entry at \a r x \a c - Entry entry(size_t r, size_t c, Entry const& e) { - self()->entries_[self->index(r, c)] = e; - return entry(r, c); - } - - //! @brief Get the entry at \a r x \a c - EntryRef entryGet(size_t r, size_t c) { - return self()->entries_[self->index(r, c)]; - } - - /*! - * @brief Change the entries from \a rFirst x \a cFirst to \a rLast x \a cLast - * - * @param [in] rFirst - * @param [in] rLast - * @param [in] cFirst - * @param [in] cLast - * @param [in] e value - * @return void - */ - void entries(ssize_t rFirst, ssize_t rLast, - ssize_t cFirst, ssize_t cLast, - Entry const& e) { - for (ssize_t r = rFirst; r <= rLast; r++) { - for (ssize_t c = cFirst; c <=cFirst; c++) { - entry(r, c, e); - } - } - } - - /*! - * @brief Return a \a rLast-rFirst+1 x \a cLast-cFirst+1 matrix - * - * With value be the entries from \a rFirst x \a cFirst to \a rLast x \a cLast - * - * @param [in] rFirst - * @param [in] rLast - * @param [in] cFirst - * @param [in] cLast - * @return a matrix - */ - Matrix subMatrix(size_t rFirst, size_t rLast, - size_t cFirst, size_t cLast) const { - if (rFirst > rLast || cFirst > cLast) return Matrix(); - if (rFirst == 0 && cFirst == 0) { - Matrix ret(*this); - ret.size(rLast + 1, cLast + 1, Entry(0)); - return ret; - } - Matrix ret(rLast - rFirst + 1, cLast - cFirst + 1, entry(rFirst, cFirst)); - for (size_t r = rFirst; r <= rLast; r++) - for (size_t c = cFirst; c <= cLast; c++) - ret.entry(r - rFirst, c - cFirst, entry(r, c)); - return ret; - } - - //! @brief Return the \a r -th row - Matrix row(size_t r) const { - return subMatrix(r, r, 0, cols() - 1); - } - - //! @brief Return the \a c -th column - Matrix col(size_t c) const { - return subMatrix(0, rows() - 1, c, c); - } - - //! @brief return +\a (*this) - Matrix positive() const { - return *this; - } - - //! @brief return -\a (*this) - Matrix negative() const { - Matrix ret(*this); - for (size_t r = 0, R = rows(); r < R; r++) - for (size_t c = 0, C = cols(); c < C; c++) - ret.entry(r, c, -ret.entry(r, c)); - return ret; - } - - /*! @brief return \a (*this) + \a m. - * - * If the size not match, it will return an invalid matrix - */ - Matrix add(Matrix const& m) const { - if (rows() != m.rows() || cols() != m.cols()) return Matrix(); - Matrix ret(*this); - for (size_t r = 0, R = rows(); r < R; r++) - for (size_t c = 0, C = cols(); c < C; c++) - ret.entry(r, c, ret.entry(r, c) + m.entry(r, c)); - return ret; - } - - /*! @brief return \a (*this) - \a m. - * - * If the size not match, it will return an invalid matrix - */ - Matrix sub(Matrix const& m) const { - if (rows() != m.rows() || cols() != m.cols()) return Matrix(); - Matrix ret(*this); - for (size_t r = 0, R = rows(); r < R; r++) - for (size_t c = 0, C = cols(); c < C; c++) - ret.entry(r, c, ret.entry(r, c) - m.entry(r, c)); - return ret; - } - - /*! @brief return \a (*this) times \a m. - * - * If the size not match, it will return an invalid matrix - */ - Matrix mul(Matrix const& m) const { - if (cols() != m.rows()) return Matrix(); - Matrix ret(rows(), m.cols(), Entry(0)); - for (size_t r = 0, R = rows(); r < R; r++) - for (size_t c = 0, C = m.cols(); c < C; c++) - for (size_t k = 0, K = cols(); k < K; k++) - ret.entry(r, c, ret.entry(r, c) + entry(r, k) * m.entry(k, c)); - return ret; - } - - //! @brief return \a (*this) times \a s. \a s is a scalar - Matrix mul(Entry const& s) const { - Matrix ret(*this); - for (size_t r = 0, R = rows(); r < R; r++) - for (size_t c = 0, C = cols(); c < C; c++) - ret.entry(r, c, ret.entry(r, c) * s); - return ret; - } - - //! @brief return \a (*this) / \a s. \a s is a scalar - Matrix div(Entry const& s) const { - Matrix ret(*this); - for (size_t r = 0, R = rows(); r < R; r++) - for (size_t c = 0, C = cols(); c < C; c++) - ret.entry(r, c, ret.entry(r, c) / s); - return ret; - } - - //! @brief Return a identity matrix with size equal to itself - Matrix identity() const { - Matrix ret(*this); - ret.identitied(); - return ret; - } - - /*! - * @brief Let itself be an identity matrix - * - * Our definition of Identity matrix is 1 for entry(i, i) and 0 otherwise. - */ - Matrix& identitied() { - for (size_t r = 0, R = rows(); r < R; r++) - for (size_t c = 0, C = cols(); c < C; c++) - entry(r, c, (r == c ? Entry(1) : Entry(0))); - return *this; - } - - /*! - * @brief Let itself be an diagonal form of original itself - */ - Matrix& diagonaled() { - triangulared(); - for (size_t i = 0, I = rows(); i < I; ++i) { - for (size_t j = i + 1, J = cols(); j < J; ++j) { - entry(i, j, Entry(0)); - } - } - return *this; - } - - /*! - * @brief Return a matrix which is a diangonal form of me - */ - Matrix diagonal() const { - Matrix ret(*this); - ret.diagonaled(); - return ret; - } - - /*! - * @brief Return a matrix which is an inverse matrix of \a (*this) - * - * If inverse matrix doesn't exist, it will return a invalid matrix - */ - Matrix inverse() const { - if (rows() != cols() || rows() == 0) return Matrix<Entry>(); - Matrix tmp(rows(), cols() * 2, Entry(0)); - for (size_t r = 0, R = rows(); r < R; r++) { - for (size_t c = 0, C = cols(); c < C; c++) { - tmp.entry(r, c, entry(r, c)); - tmp.entry(r, c + cols(), (r == c ? Entry(1) : Entry(0))); - } - } - tmp.triangulared(); - for (ssize_t r = rows() - 1; r >= 0; r--) { - if (tmp(r, r) == Entry(0)) return Matrix<Entry>(); - for (ssize_t r2 = r - 1; r2 >= 0; r2--) { - Entry rat(-tmp.entry(r2, r) / tmp.entry(r, r)); - for (size_t c = r, C = tmp.cols(); c < C; c++) { - tmp.entry(r2, c, tmp.entry(r2, c) + rat * tmp(r, c)); - } - } - Entry rat(tmp.entry(r, r)); - for (size_t c = cols(), C = tmp.cols(); c < C; c++) { - tmp.entry(r, c - cols(), tmp.entry(r, c) / rat); - } - } - tmp.size(cols(), rows(), Entry(0)); - return tmp; - } - - //! @brief let itself become itself's inverse matrix - Matrix& inversed() { - copyFrom(inverse()); - return *this; - } - - //! @brief return itself's transpose matrix - Matrix transpose() const { - Matrix ret(cols(), rows(), Entry(0)); - for (size_t r = 0, R = cols(); r < R; r++) - for (size_t c = 0, C = rows(); c < C; c++) - ret.entry(r, c, entry(c, r)); - return ret; - } - - //! @brief Let itself become itself's transpose matrix - Matrix& transposed() { - copyFrom(transpose()); - return *this; - } - - //! @brief return a matrix which is the triangular form of \a (*this) - Matrix triangular() const { - Matrix<Entry> ret(*this); - ret.triangulared(); - return ret; - } - - //! @brief triangluar itself - Matrix& triangulared() { - for (size_t r = 0, c = 0, R = rows(), C = cols(); r < R && c < C; r++) { - ssize_t maxR; - for ( ; c < C; c++) { - maxR = -1; - for (size_t r2 = r; r2 < R; r2++) - if (maxR == -1 || tAbs(entry(r2, c)) > tAbs(entry(maxR, c))) - maxR = r2; - if (entry(maxR, c) != Entry(0)) break; - } - if (c >= C) break; - if (maxR != (ssize_t)r) { - for (size_t c2 = c; c2 < C; c2++) - std::swap(self()->entries_[self->index( r, c2)], - self()->entries_[self->index(maxR, c2)]); - } - for (size_t r2 = r + 1; r2 < R; r2++) { - Entry rati = -entry(r2, c) / entry(r, c); - entry(r2, c, Entry(0)); - for (size_t c2 = c + 1; c2 < C; c2++) - entry(r2, c2, entry(r2, c2) + entry(r, c2) * rati); - } - } - return *this; - } - - //! @brief same as \a copyFrom - Matrix& operator=(Matrix const& m) { - return copyFrom(m); - } - - //! @brief same as \a entry(r,c) - Entry operator()(size_t r, size_t c) const { - return entry(r, c); - } - - //! @brief same as \a entry(r,c,e) - Entry operator()(size_t r, size_t c, Entry const& e) { - return entry(r, c, e); - } - - //! @brief same as \a positive() - Matrix operator+() const { - return positive(); - } - - //! @brief same as \a negative() - Matrix operator-() const { - return negative(); - } - - //! @brief same as \a add(m) - Matrix operator+(Matrix const& m) const { - return add(m); - } - - //! @brief same as \a sub(m) - Matrix operator-(Matrix const& m) const { - return sub(m); - } - - //! @brief same as \a mul(m) - Matrix operator*(Matrix const& m) const { - return mul(m); - } - - //! @brief same as \a mul(m) - Matrix operator*(Entry const& s) const { - return mul(s); - } - - //! @brief same as \a div(s) - Matrix operator/(Entry const& s) const { - return div(s); - } -}; - -} // meow - -#endif // math_Matrix_H__ diff --git a/meowpp/math/Transformation.h b/meowpp/math/Transformation.h deleted file mode 100644 index 086cb03..0000000 --- a/meowpp/math/Transformation.h +++ /dev/null @@ -1,237 +0,0 @@ -#ifndef math_Transformation_H__ -#define math_Transformation_H__ - -#include "Matrix.h" -#include "../Self.h" - -#include <list> -#include <cstdlib> - -namespace meow { - -/*! - * @brief A base class for implementing kinds of transformations. - * - * We define that the input and output form of our transformations all be - * \b matrix . Some advance methods such as calculating jacobian matrix - * will require that the input form must be a vector. - * @author cat_leopard - */ -template<class Scalar> -class Transformation { -private: - struct Myself { - size_t inputRows_; - size_t inputCols_; - size_t outputRows_; - size_t outputCols_; - size_t psize_; - - Myself(Myself const& b): - inputRows_(b.inputRows_), inputCols_(b.inputCols_), - outputRows_(b.outputRows_), outputCols_(b.outputCols_), - psize_(b.psize_) { - } - - Myself(size_t ir, size_t ic, size_t or_, size_t oc, size_t ps): - inputRows_(ir), inputCols_(ic), outputRows_(or_), outputCols_(oc), - psize_(ps) { - } - - ~Myself() { - } - }; - - Self<Myself> const self; -protected: - /*! - * Construct and setup - * @param [in] inputRows number of rows of the input matrix. - * @param [in] inputCols number of columns of the input matrix. - * @param [in] outputRows number of rows of the output matrix. - * @param [in] outputCols number of columns of the output matrix. - * @param [in] psize number of parameters - */ - Transformation(size_t inputRows, size_t inputCols, - size_t outputRows, size_t outputCols, - size_t psize): - self(Myself(inputRows, inputCols, outputRows, outputCols, psize)) { - } - - /*! - * Construct and copy setings from another transformation class. - * @param [in] b Specify where to copy the informations. - */ - Transformation(Transformation const& b): - self(b.self, Self<Myself>::COPY_FROM) { - } - - /*! - * @brief Copy from the specified one - * - * @param [in] b The specified one - * @return \c *this - */ - Transformation& copyFrom(Transformation const& b) { - self().copyFrom(b.self); - return *this; - } - - /*! - * @brief reference from the specified one - * - * @param [in] b The specified one - * @return \c *this - */ - Transformation& referenceFrom(Transformation const& b) { - self().referenceFrom(b.self); - return *this; - } -public: - /*! - * Destructor - */ - virtual ~Transformation() { - } - - /*! - * @brief Return the number of rows of the input matrix. - * - * @return Number of rows. - */ - size_t inputRows() const { - return self->inputRows_; - } - - /*! - * @brief Return the number of columns of the input matrix. - * - * @return Number of columns. - */ - size_t inputCols() const { - return self->inputCols_; - } - - /*! - * @brief Return the number of rows of the output matrix. - * - * @return Number of rows. - */ - size_t outputRows() const { - return self->outputRows_; - } - - /*! - * @brief Return the number of columns of the output matrix. - * - * @return Number of columns. - */ - size_t outputCols() const { - return self->outputCols_; - } - - /*! - * @brief Return the number of parameters. - * - * @return Number of parameters. - */ - size_t parameterSize() const { - return self->psize_; - } - - /*! - * @brief Get the \a i -th parameter. - * - * @param [in] i The index of the specified parameter. - * @note It's a pure virtual method. - */ - virtual Scalar parameter(size_t i) const = 0; - - /*! - * @brief Setup the \a i -th parameter. - * - * @param [in] i The index of the specified parameter. - * @param [in] s The new value to the specified parameter. - * @note It's a pure virtual method. - */ - virtual Scalar parameter(size_t i, Scalar const& s) = 0; - - /*! - * @brief Do transformate. - * - * @param [in] x The input matrix. - * @note It's a pure virtual method. - */ - virtual Matrix<Scalar> transformate(Matrix<Scalar> const& x) const = 0; - - /*! - * @brief Calculate the jacobian matrix (derivate by the input matrix) - * of the transformation. - * - * Consider the case of a non-differentiable - * transformation might be implemented, we return an empty matrix - * now instead of making it be a pure virtual method. - * @param [in] x The input matrix. - * @return An empty matrix. - */ - virtual Matrix<Scalar> jacobian(Matrix<Scalar> const& x) const { - return Matrix<Scalar>(); - } - - /*! - * @brief Calculate the jacobian matrix (derivate by the \a i -th parameter) - * of the transformation. - * - * Consider the case of a non-differentiable transformation might be - * implemented, we return an empty matrix now instead of making it be - * a pure virtual method. - * @param [in] x The input matrix. - * @param [in] i The index of the specified parameter. - * @return An empty matrix. - */ - virtual Matrix<Scalar> jacobian(Matrix<Scalar> const& x, size_t i) const { - return Matrix<Scalar>(); - } - - /*! - * @brief Return whether this transformation is inversable or not - * - * @return \c false - */ - virtual bool inversable() const { return false; } - - /*! - * @brief Do the inverse transformation - * - * @param [in] x The input matirx - * @return An empty matrix - */ - virtual Matrix<Scalar> transformateInv(Matrix<Scalar> const& x) const { - return Matrix<Scalar>(); - } - - /*! - * @brief Return the jacobian matrix of the inverse transformation - * - * @param [in] x The input matirx - * @return An empty matrix - */ - virtual Matrix<Scalar> jacobianInv(Matrix<Scalar> const& x) const { - return Matrix<Scalar>(); - } - - /*! - * @brief Return the jacobian matrix of the inverse transformation - * - * @param [in] x The input matirx - * @param [in] i The index of the specified parameter. - * @return An empty matrix - */ - virtual Matrix<Scalar> jacobianInv(Matrix<Scalar> const& x, size_t i) const { - return Matrix<Scalar>(); - } -}; - -} // meow - -#endif // math_Transformation_H__ diff --git a/meowpp/math/Transformations.h b/meowpp/math/Transformations.h deleted file mode 100644 index 99d6483..0000000 --- a/meowpp/math/Transformations.h +++ /dev/null @@ -1,550 +0,0 @@ -#ifndef math_Transformations_H__ -#define math_Transformations_H__ - -#include "Transformation.h" -#include "Matrix.h" -#include "utility.h" -#include "../Self.h" - -#include <cstdlib> - -namespace meow { - -/*! - * @brief A ball projection is to project the given vector to a hyper-sphere - * - * Assume: - * - The dimension of a ball projection is \f$ N \f$ - * - The radius of the hyper-sphere is \f$ R \f$ - * . - * Then the transformation is like below: \n - * \f[ - * \left[ - * \begin{array}{c} - * x_1 \\ - * x_2 \\ - * x_3 \\ - * . \\ - * . \\ - * . \\ - * x_N \\ - * \end{array} - * \right] - * \stackrel{transformate}{\rightarrow} - * \left[ - * \begin{array}{c} - * \frac{x_1 \times R}{L} \\ - * \frac{x_2 \times R}{L} \\ - * \frac{x_3 \times R}{L} \\ - * . \\ - * . \\ - * . \\ - * \frac{x_N \times R}{L} \\ - * \end{array} - * \right] \\ - * \f] - * where \f$ L=\sqrt{x_1^2 + x_2^2 + x_3^2 + ... + x_N^2 } \f$ - * @author cat_leopard - */ -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: - /*! - * Constructor, copy settings from given BallProjection - * @param [in] b another ball projection class - */ - BallProjection(BallProjection const& b): Transformation<Scalar>(b), - self(b.self, Self<Myself>::COPY_FROM) { - } - - /*! - * Constructor and setup, radius = 1 - * @param [in] d Dimension of the input/output vector - */ - BallProjection(size_t d): Transformation<Scalar>(d, 1, d, 1, 1), - self(Myself(d)) { - radius(1); - } - - /*! - * Constructor and setup - * @param [in] d Dimension of the input/output vector - * @param [in] r Radius of the hyper-sphere - */ - BallProjection(size_t d, Scalar const& r): Transformation<Scalar>(d,1,d,1,1), - self(Myself(d, r)) { - radius(r); - } - - /*! - * @brief Copy settings from another one - * @param [in] b Another one - * @return \c *this - */ - BallProjection& copyFrom(BallProjection const& b) { - Transformation<Scalar>::copyFrom(b); - copyFrom(b); - return *this; - } - - /*! - * @brief Reference settings from another one - * @param [in] b Another one - * @return \c *this - */ - BallProjection& referenceFrom(BallProjection const& b) { - Transformation<Scalar>::referenceFrom(b); - referenceFrom(b); - return *this; - } - - /*! - * @brief same as \c radius() - */ - Scalar parameter(size_t i) const { - return radius(); - } - - /*! - * @brief same as \c radius(s) - */ - Scalar parameter(size_t i, Scalar const& s) { - return radius(s); - } - - /*! - * @brief Return the value of the radius - */ - Scalar radius() const { - return self->radius_; - } - - /*! - * @brief Setup the radius - * - * @param [in] r New value of the radius - * @return New radius - */ - Scalar radius(Scalar const& r) { - self()->radius_ = r; - return radius(); - } - - /*! - * @brief Get the dimension of this projection - */ - size_t dimension() const { - return self->dimension_; - } - - - /*! - * @brief Project the input vector(s) onto the hyper-sphere and return it. - * - * If the number of columns of the input matrix is larger than 1, this - * method will think that you want to transform multiple vector once - * and the number of columns of the output matrix will be the same of - * the number of columns of the input one. - * - * @param [in] x The input matrix. - * @return The output matrix. - * @note Take into account that too much safty checking will lead to - * inefficient, this method will not checking whether the dimension - * of the input vector/matrix is right. So be sure the data is valid - * before you call this method. - */ - 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; - } - - /*! - * @brief Return the jacobian matrix (derivate by the input vector) - * of this projection. - * - * This method only allow a vector-like matrix be input. - * Assume: - * - The dimension of a ball projection is \f$ N \f$ - * - The length of the input vector is \f$ L=\sqrt{x_1^2+x_2^2+...+x_N^2} \f$ - * - The radius of the hyper-sphere is \f$ R \f$ - * . - * Then the jacobian matrix is like below: \n - * \f[ - * \frac{R}{L^3} \times \left[ - * \begin{array}{ccccc} - * L^2-x_1^2 & -x_1x_2 & -x_1x_3 & ... & -x_1x_N \\ - * -x_2x_1 & L^2-x_2^2 & -x_2x_3 & ... & -x_2x_N \\ - * -x_3x_1 & -x_3x_2 & L^2-x_3^2 & ... & -x_3x_N \\ - * . & . & . & & . \\ - * . & . & . & & . \\ - * . & . & . & & . \\ - * -x_Nx_1 & -x_Nx_2 & -x_Nx_3 & ... & L^2-x_N^2 \\ - * \end{array} - * \right] - * \f] - * - * @param [in] x The input matrix. - * @return The output matrix. - */ - 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; - } - - /*! - * @brief Return the jacobian matrix (derivate by radius) of this projection. - * - * This method only allow a vector-like matrix be input. - * Assume: - * - The dimension of a ball projection is \f$ N \f$ - * - The length of the input vector is \f$ L=\sqrt{x_1^2+x_2^2+...+x_N^2} \f$ - * - The radius of the hyper-sphere is \f$ R \f$ - * . - * Then the jacobian matrix is like below: \n - * \f[ - * R \times \left[ - * \begin{array}{c} - * \frac{x_1}{L} \\ - * \frac{x_2}{L} \\ - * \frac{x_3}{L} \\ - * . \\ - * . \\ - * . \\ - * \frac{x_N}{L} \\ - * \end{array} - * \right] - * \f] - * - * @param [in] x The input matrix. - * @param [in] i Useless parameter - * @return The output matrix. - */ - 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))); - } - - /*! - * @brief Same as \c copyFrom(b) - */ - BallProjection& operator=(BallProjection const& b) { - return copyFrom(b); - } - - /*! - * @brief Same as \c transformate(v) - */ - Matrix<Scalar> operator()(Matrix<Scalar> const& v) const { - return transformate(v); - } -}; - - -/*! - * @brief A \b photo \b projection is a kind of transformation that project - * point/vector to a flat \b photo - * - * Assume: - * - The dimension of a photo projection is \f$ N \f$ - * - The length of the input vector is \f$ L \f$ - * - The focal length is \f$ f \f$ - * . - * Then transformation is like below: \n - * \f[ - * \left[ - * \begin{array}{c} - * x_1 \\ - * x_2 \\ - * x_3 \\ - * . \\ - * . \\ - * . \\ - * x_N \\ - * \end{array} - * \right] - * \stackrel{transformate}{\rightarrow} - * \left[ - * \begin{array}{c} - * \frac{-x_1 \times f}{x_N} \\ - * \frac{-x_2 \times f}{x_N} \\ - * \frac{-x_3 \times f}{x_N} \\ - * . \\ - * . \\ - * . \\ - * -f \\ - * \end{array} - * \right] \\ - * \f] - * i.e. projecte the vector onto the plane \f$ x_N = -f \f$. - * - * @author cat_leopard - */ -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: - /*! - * Constructor, focal = 1 - */ - PhotoProjection(size_t dimension): - Transformation<Scalar>(dimension, 1, dimension, 1, 1), - self(Myself(dimension, 1)) { - } - - /*! - * Constructor - */ - PhotoProjection(size_t dimension, Scalar const& f): - Transformation<Scalar>(dimension, 1, dimension, 1, 1), - self(Myself(dimension, f)) { - } - - /*! - * Constructor, copy settings from another PhotoProjection. - */ - PhotoProjection(PhotoProjection const& p): Transformation<Scalar>(p), - self(p.self, Self<Myself>::COPY_FROM) { - } - - /*! - * Copy settings from another one - * @param [in] b another one - * @return \c *this - */ - PhotoProjection& copyFrom(PhotoProjection const& b) { - Transformation<Scalar>::copyFrom(b); - self().copyFrom(b.self); - return *this; - } - - /*! - * Reference settings from another one - * @param [in] b another one - * @return \c *this - */ - PhotoProjection& referenceFrom(PhotoProjection const& b) { - Transformation<Scalar>::referenceFrom(b); - self().referenceFrom(b.self); - return *this; - } - - /*! - * @brief Same as \c focal() - */ - Scalar parameter(size_t i) const { - return focal(); - } - - /*! - * @brief Same as \c focal(s) - */ - Scalar parameter(size_t i, Scalar const& s){ - return focal(s); - } - - /*! - * @brief Get the focal length - * @return Focal length - */ - Scalar focal() const { - return self->focal_; - } - - /*! - * @brief Set the focal length - * - * @param [in] f New focal length - * @return New focal length - */ - Scalar focal(Scalar const& f){ - self()->focal_ = f; - return focal(); - } - - /*! - * @brief Get the dimension of this projection - */ - size_t dimension() const { - return self->dimension_; - } - - /*! - * @brief Project the input vector(s) onto the plane - * - * The equation of the plane is \f$ x_N = -f \f$, where the \f$ N \f$ - * is the dimension of this projection and f is the focal length. \n - * If the number of columns of the input matrix is larger than 1, this - * method will think that you want to transform multiple vector once - * and the number of columns of the output matrix will be the same of - * the number of columns of the input one. - * - * @param [in] x The input matrix. - * @return The output matrix. - * @note Take into account that too much safty checking will lead to - * inefficient, this method will not checking whether the dimension - * of the input vector/matrix is right. So be sure the data is valid - * before you call this method. - */ - 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; - } - - /*! - * @brief Return the jacobian matrix (derivate by the input vector) - * of this projection. - * - * This method only allow a vector-like matrix be input. - * Assume: - * - The dimension of this projection is \f$ N \f$ - * - The length of the input vector is \f$ L=\sqrt{x_1^2+x_2^2+...+x_N^2} \f$ - * - The focal length of this projection is \f$ f \f$ - * . - * Then the jacobian matrix is like below: \n - * \f[ - * f \times - * \left[ - * \begin{array}{ccccc} - * \frac{-1}{x_N} & 0 & 0 & ... & \frac{1}{x_N^2} \\ - * 0 & \frac{-1}{x_N} & 0 & ... & \frac{1}{x_N^2} \\ - * 0 & 0 & \frac{-1}{x_N} & ... & \frac{1}{x_N^2} \\ - * . & . & . & & . \\ - * . & . & . & & . \\ - * . & . & . & & . \\ - * 0 & 0 & 0 & ... & 0 \\ - * \end{array} - * \right] - * \f] - * - * @param [in] x The input matrix. - * @return The output matrix. - */ - 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; - } - - /*! - * @brief Return the jacobian matrix (derivate by the focus length) - * of this projection. - * - * This method only allow a vector-like matrix be input. - * Assume: - * - The dimension of this projection is \f$ N \f$ - * - The length of the input vector is \f$ L=\sqrt{x_1^2+x_2^2+...+x_N^2} \f$ - * - The focal length of this projection is \f$ f \f$ - * . - * Then the jacobian matrix is like below: \n - * \f[ - * \left[ - * \begin{array}{c} - * \frac{-x_1}{x_N} \\ - * \frac{-x_2}{x_N} \\ - * \frac{-x_3}{x_N} \\ - * . \\ - * . \\ - * . \\ - * -1 \\ - * \end{array} - * \right] - * \f] - * - * @param [in] x The input matrix. - * @param [in] i Useless parameter - * @return The output matrix. - */ - 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; - } - - /*! - * @brief Same as \c copyFrom(b) - */ - PhotoProjection& operator=(PhotoProjection const& b) { - return copyFrom(b); - } - - /*! - * @brief Same as \c transformate(v) - */ - Matrix<Scalar> operator()(Matrix<Scalar> const& v) const { - return transformate(v); - } -}; - -} // meow - -#endif // Transformations_H__ diff --git a/meowpp/math/Vector.h b/meowpp/math/Vector.h deleted file mode 100644 index f72b043..0000000 --- a/meowpp/math/Vector.h +++ /dev/null @@ -1,267 +0,0 @@ -#ifndef math_Vector_H__ -#define math_Vector_H__ - -#include "../Self.h" -#include "Matrix.h" - -#include <vector> - -#include <cmath> - -namespace meow { - -/*! - * @brief \b vector - * - * @author cat_leopard - */ -template<class Scalar> -class Vector { -public: - typedef typename Matrix<Scalar>::EntryRefK ScalarRefK; - typedef typename Matrix<Scalar>::EntryRef ScalarRef ; -private: - Matrix<Scalar> matrix_; -public: - /*! - * @brief constructor - * - * With \b dimension=0, which means \b invalid. - */ - Vector() { - } - - /*! - * @brief constructor - * - * Copy from another vector - * - * @param [in] v another vector - */ - Vector(Vector const& v): matrix_(v.matrix_) { - } - - /*! - * @brief constructor - * - * From matrix's first column - * - * @param [in] m matrix - */ - Vector(Matrix<Scalar> const& m): matrix_(m.col(0)) { - } - - /*! - * @brief constructor - * - * Copy from another std::vector - * - * @param [in] v vector - */ - Vector(std::vector<Scalar> const& v): matrix_(v.size(), 1, Scalar(0)) { - for (size_t i = 0, I = v.size(); i < I; i++) { - matrix_.entry(i, 0, v[i]); - } - } - - /*! - * @brief constructor - * - * setup dimension and inital value - * - * @param [in] d dimension - * @param [in] e inital value - */ - Vector(size_t d, Scalar const& e): matrix_(d, 1, e) { - } - - //! @brief destructor - ~Vector() { - } - - //! @brief copy from ... - Vector& copyFrom(Vector const& v) { - matrix_.copyFrom(v.matrix_); - return *this; - } - - //! @brief reference from ... - Vector& referenceFrom(Vector const& v) { - matrix_.referenceFrom(v.matrix_); - return *this; - } - - //! @brief Return a \a dimension x 1 matrix form of it - Matrix<Scalar> matrix() const { - return matrix_; - } - - //! @brief return dimension - size_t dimension() const { - return matrix_.rows(); - } - - /*! - * @brief resize the dimension - * - * @param [in] d new dimension - * @param [in] s inital entry - * @return new dimension - */ - size_t dimension(size_t d, Scalar const& s) { - matrix_.rows(d, s); - return dimension(); - } - - /*! - * @brief Return whether \c dimension>0 is true or not - * @return \c true/false - */ - bool valid() const { - return (dimension() > 0); - } - - //! @brief return \a i -th scalar - Scalar scalar(size_t i) const { - return matrix_.entry(i, 0); - } - - /*! - * @brief change \a i -th scalar - * - * @param [in] i i-th - * @param [in] s new value - */ - Scalar scalar(size_t i, Scalar const& s) { - matrix_.entry(i, 0, s); - return scalar(i); - } - - //! @brief return \a i -th scalar with non-constant type - ScalarRef scalarGet(size_t i) { - return matrix_.entryGet(i); - } - - /*! - * @brief change \a i -th to \a j -th scalars - * - * @param [in] i i-th - * @param [in] j j-th - * @param [in] s new value - */ - void scalars(size_t i, size_t j, Scalar const& s) { - for (size_t it = i; it <= j; ++it) { - matrix_.entry(it, 0, s); - } - } - - //! @brief subvector form i-th to j-th - Vector subVector(size_t i, size_t j) { - return Vector(matrix_.subMatrix(i, 0, j, 0)); - } - - //! @brief return +\a (*this) - Vector positive() const { - return *this; - } - - //! @brief return -\a (*this) - Vector negative() const { - return Vector(matrix_.negative()); - } - - //! @brief return \a (*this)+v - Vector add(Vector const& v) const { - return Vector(matrix_.add(v.matrix_)); - } - - //! @brief return \a (*this)-v - Vector sub(Vector const& v) const { - return Vector(matrix_.sub(v.matrix_)); - } - - //! @brief return \a (*this)*s , where s is a scalar - Vector mul(Scalar const& s) const { - return Vector(matrix_.mul(s)); - } - - //! @brief return \a (*this)/s , where s is a scalar - Vector div(Scalar const& s) const { - return Vector(matrix_.div(s)); - } - - //! @brief dot - Scalar dot(Vector const& v) const { - return matrix_.transpose().mul(v.matrix_).entry(0, 0); - } - - //! @brief sqrt of \a length2 - Scalar length() const { - return Scalar(sqrt((double)length2())); - } - - //! @brief same as \a (*this).dot(*this) - Scalar length2() const { - return dot(*this); - } - - //! @brief return a normalize form of itself - Vector normalize() const { - return div(length()); - } - - //! @brief Let itself be normalize form - Vector& normalized() { - copyFrom(normalize()); - return *this; - } - - //! @brief same as copyFrom - Vector& operator=(Vector const& v) { - return copyFrom(v); - } - - //! @brief same as entry(i) - Scalar operator()(size_t i) const { - return scalar(i); - } - - //! @brief same as positive() - Vector operator+() const { - return positive(); - } - - //! @brief same as negative() - Vector operator-() const { - return negative(); - } - - //! @brief same as add(v) - Vector operator+(Vector const& v) const { - return add(v); - } - - //! @brief same as sub(v) - Vector operator-(Vector const& v) const { - return sub(v); - } - - //! @brief same as dot(v) - Scalar operator*(Vector const& v) const { - return dot(v); - } - - //! @brief same as mul(s) - Vector operator*(Scalar const& s) const { - return mul(s); - } - - //! @brief same as div(s) - Vector operator/(Scalar const& s) const { - return div(s); - } -}; - -} // meow - -#endif // math_Vector_H__ diff --git a/meowpp/math/methods.h b/meowpp/math/methods.h deleted file mode 100644 index 06125b5..0000000 --- a/meowpp/math/methods.h +++ /dev/null @@ -1,230 +0,0 @@ -#ifndef math_methods_H__ -#define math_methods_H__ - -#include "Matrix.h" -#include "Vector.h" -#include "utility.h" - -#include <cstdlib> -#include <vector> - -namespace meow { - -/*! - * @brief Run the \b RANSAC method to approach the best solution. - * - * \b RANdom \b SAmple \b Consensus is an iterative method to estimate - * parameters of a mathematical model from a set of observed data which - * contains \c outliers. \n - * Each iterator it will choose a subset of elements, the smallest set which can - * form a valid parameters, from the data set. And then calculate how many - * elements in the whole data set is inliers. After iterator much times, - * we just say the best solution is the parameters that has the much - * inliers elements in whole iterators. - * - * Assume: - * - We need at least \f$ N \f$ element to form a valid parameters. - * - The probability of choosing a right element from data set each time is - * \f$ p_0 \f$. - * - We want the probability of our solution actually being the best solution - * be \f$ P \f$. - * - We need to iterator \f$ M \f$ times. - * . - * Then we can estimate the number of iterations \f$ M \f$ : - * \f[ - * \begin{aligned} - * & (1 - p_0^N)^M \leq(1 - P) \\ - * \Rightarrow & M \log(1 - p_0^N) \leq \log(1 - P) \\ - * \Rightarrow & M \geq \frac{\log(1 - p)}{\log(1 - p_0^N)},~~ - * \because (1-p_0^N<1 \Rightarrow \log(1-p_0^N)<0) - * \end{aligned} - * \f] - * - * So in this function we choose - * \f$ M = \lceil \frac{\log(1 - P)}{\log(1 - p_0^N)} \rceil \f$ - * - * @param [in] data The whole data sett - * @param [in] w Weight function to give a floating number for a given - * parameters which means how best this solution is. Negitave - * number means invalid parameters. - * @param [in] N \f$ N \f$, defined above - * @param [in] p0 \f$ p_0 \f$, defined above - * @param [in] P \f$ P \f$, defined above - * @return solution. - * - * @author cat_leopard - */ -template<class Data, class WeightingClass> -inline std::vector<Data> ransac(std::vector<Data> const& data, - WeightingClass const& w, - size_t N, - double p0, double P) { - if (data.size() < N) { - return std::vector<Data>(); - } - double ww = -1.0; - std::vector<Data> ret; - for (double count = ceil(log(1.0 - P) / log(1.0 - pow(p0, N))); - count > 0.0; count -= 1.0) { - std::vector<Data> sample; - std::vector<int> index(N); - for (size_t i = 0; i < N; i++) { - for (bool ok = false; !ok; ) { - index[i] = rand() % data.size(); - ok = true; - for (size_t j = 0; ok && j < i; j++) - if (index[i] == index[j]) - ok = false; - } - sample.push_back(data[index[i]]); - } - double w_now = w(sample, data); - if (w_now < 0) { - count += 0.5; - continue; - } - if (ww < w_now) { - ret = sample; - ww = w_now; - } - } - return ret; -} - - -/* - * @brief Run the \b Levenberg-Marquardt method to solve a non-linear - * least squares problem. - * - * Assume: - * - The function we want to optimize is - * \f$ F: \mathbb{R} ^N \mapsto \mathbb{R}^M \f$ - * - We want to find the best solution \f$ v \f$ such that - * \f$ F(v)^T F(v) = 0\f$. But there is a gived threshold - * \f$ \epsilon \f$, we can just find a \f$ v \f$ such that - * \f$ F(v)^T F(v) < \epsilon \f$, which is mush easier. - * - User gived a initiial vector \f$ v_0 \f$ - * . - * Then we just iteratilly find \f$ v_1, v_2, v_3, v_4... \f$ until a - * vector \f$ v_k \f$ satisified that \f$ F(v_k)^TF(v_k)<\epsilon \f$ . - * And each iterator we have: - * \f[ - * v_{i+1} = v_i + (J(v_i)^TJ(v_i)+\lambda I_{N\times N})^{-1} J(v_i)^T F(v_i) - * \f] - * Where \f$ J(v) \f$ is a jacobian matrix defined below: - * \f[ - * J(v) = \frac{d}{dv}F(v) = - * \left[ \begin{array}{ccccc} - * \frac{\partial F_1(v)}{\partial v_1} & - * \frac{\partial F_1(v)}{\partial v_2} & - * \frac{\partial F_1(v)}{\partial v_3} & - * ... & - * \frac{\partial F_1(v)}{\partial v_N} \\ - * \frac{\partial F_2(v)}{\partial v_1} & - * \frac{\partial F_2(v)}{\partial v_2} & - * \frac{\partial F_2(v)}{\partial v_3} & - * ... & - * \frac{\partial F_2(v)}{\partial v_N} \\ - * \frac{\partial F_3(v)}{\partial v_1} & - * \frac{\partial F_3(v)}{\partial v_2} & - * \frac{\partial F_3(v)}{\partial v_3} & - * ... & - * \frac{\partial F_3(v)}{\partial v_N} \\ - * . & . & . & & . \\ - * . & . & . & & . \\ - * . & . & . & & . \\ - * \frac{\partial F_M(v)}{\partial v_1} & - * \frac{\partial F_M(v)}{\partial v_2} & - * \frac{\partial F_M(v)}{\partial v_3} & - * ... & - * \frac{\partial F_M(v)}{\partial v_N} \\ - * \end{array} \right] - * \f] - * And \f$ \lambda \f$ is a magic number.... - * @param [in] func \f$ F \f$, a function(class with \c operator() ) - * which input a vector and the output the squares errors. - * @param [in] jaco \f$ J \f$, a function which input a vector - * and then output \b func derivate by the vector - * @param [in] iden \f$ \lambda I_{N \times N} \f$, defined above - * @param [in] init \f$ v_0 \f$Initial vector - * @param [in] stop A function return a boolean which means the error is - * acceptable or not, so \f[ - * S_{top}(v) = \begin{cases} - * true & if~F(v)<\epsilon \\ - * false & else - * \end{cases} - * \f] - * @param [in] counter To prevent infinit loop. - * @return a vector which means the best solution this function found. - * - * @author cat_leopard - */ -template<class Scalar, class Function> -inline Vector<Scalar> levenbergMarquardt(Function const& f, - Vector<Scalar> const& init, - int counter = -1) { - Vector<Scalar> ans(init), residure_v; - for ( ; counter != 0 && !f.accept(residure_v = f.residure(ans)); --counter) { - Matrix<Scalar> m_j (f.jacobian(ans)); - Matrix<Scalar> m_jt(m_j.transpose()); - Matrix<Scalar> m(m_j * m_jt), M; - for (int i = 1; M.valid() == false; i++) { - M = (m + f.diagonal(ans, i)).inverse(); - } - ans = ans - M * m_jt * residure_v; - } - return ans; -} - -// residure -// jacobian -// identity -template<class Scalar, class Function> -inline Vector<Scalar> levenbergMarquardtTraining(Function & f, - Vector<Scalar> const& init, - Scalar const& init_mu, - Scalar const& mu_pow, - Scalar const& er_max, - int retry_number, - int counter) { - if (retry_number == 0) retry_number = 1; - Vector<Scalar> ans_now(init), rv_now(f.residure(ans_now)); - Vector<Scalar> ans_nxt , rv_nxt; - Scalar er_now(rv_now.length2()); - Scalar er_nxt; - Vector<Scalar> ans_best(ans_now); - Scalar er_best ( er_now); - Matrix<Scalar> m_ja, m_jt, m, iden(f.identity()); - Scalar mu(init_mu); - for ( ; counter != 0 && er_now > er_max; --counter) { - m_ja = f.jacobian(); - m_jt = m_ja.transpose(); - m = m_jt * m_ja; - bool good = false; - for (int i = 0; i != retry_number; ++i, mu = mu * mu_pow) { - ans_nxt = ans_now + (m + iden * mu).inverse() * m_jt * rv_now.matrix(); - rv_nxt = f.residure(ans_nxt); - er_nxt = rv_nxt.length2(); - if (er_nxt <= er_now) { - good = true; - break; - } - } - if (good) { - mu = mu / mu_pow; - } - mu = inRange(0.0000001, 100.0, mu); - ans_now = ans_nxt; - rv_now = rv_nxt; - er_now = er_nxt; - if (er_now < er_best) { - ans_best = ans_now; - er_best = er_now; - } - } - return ans_best; -} - -} // meow - -#endif // math_methods_H__ diff --git a/meowpp/math/utility.h b/meowpp/math/utility.h deleted file mode 100644 index 73efabb..0000000 --- a/meowpp/math/utility.h +++ /dev/null @@ -1,157 +0,0 @@ -#ifndef math_utility_H__ -#define math_utility_H__ - -#include <cstdlib> -#include <vector> -#include <algorithm> -#include <cmath> - -namespace meow { - -//! 圓周率... -static const double PI = 3.14159265358979323846264338327950288; - -/*! - * @brief 將角度調整於0~2PI - */ -template<class T> -inline T circle(T x) { - while (x < 0) x += 2.0 * PI; - while (2.0 * PI <= x) x -= 2.0 * PI; - return x; -} - -/*! - * @brief 如果abs(輸入的數值) < eps, 則回傳0, 否則回傳輸入的數值 - */ -template<class T> -inline T noEPS(T value, T eps = 1e-9) { - T epsp((eps < T(0)) ? -eps : eps); - return ((value < -epsp || value > epsp) ? value : T(0)); -} - -/*! - * @brief \c (value-lower)/(upper-lower) - */ -template<class T> -inline T normalize(T lower, T upper, T value) { - return (value - lower) / (upper - lower); -} - -/*! - * @brief \c (lower+_ratio*(upper-lower)) - */ -template<class T> -inline T denormalize(T lower, T upper, T _ratio) { - return lower + _ratio * (upper - lower); -} - -/*! - * @brief \c denormalize(l2,u2,normalize(l1,u1,m1)) - */ -template<class T> -inline T ratioMapping(T l1, T u1, T m1, T l2, T u2) { - return denormalize(l2, u2, normalize(l1, u1, m1)); -} - -/*! - * @brief \c std::min(mx,std::max(mn,v)) - */ -template<class T> -inline T inRange(T const& mn, T const& mx, T const& v) { - return std::min(mx, std::max(mn, v)); -} - -/*! - * @brief (mn <= x && x <= mx) - */ -template<class T> -inline T isInRange(T const& mn, T const& mx, T const& x) { - return (mn <= x && x <= mx); -} - -/*! - * @brief \c x*x - */ -template<class T> -inline T squ(T const& x) { - return x * x; -} - -/*! - * @brief \c x*x*x - */ -template<class T> -inline T cub(T const& x) { - return x * x * x; -} - -/*! - * @brief 只將 \c sigs 個標準差以內的數據拿來取平均 - */ -template<class T> -inline double average(T const& beg, T const& end, double sigs) { - int N = 0; - double av = 0; - for (T it = beg; it != end; ++it, ++N) { - av += *it; - } - av /= N; - double sig = 0; - for (T it = beg; it != end; ++it) { - sig += (*it - av) * (*it - av); - } - sig = sqrt(sig / N); - double lower = av - sig * sigs, upper = av + sig * sigs; - double ret = 0, retn = 0; - for (T it = beg; it != end; ++it) { - if (lower <= *it && *it <= upper) { - ret += *it; - retn++; - } - } - return ret / retn; -} - -/*! - * @brief 只將 \c sigs 個標準差以內的數據拿來取平均, 不過這次用 \c p 來加權平均 - */ -template<class T> -inline double average(T const& beg, T const& end, T const& p, double sigs) { - int N = 0; - double ps = 0; - for (T it = beg, ip = p; it != end; ++it, ++N, ++ip) { - ps += *ip; - } - double av = 0; - for (T it = beg, ip = p; it != end; ++it, ++ip) { - av += *it * *ip / ps; - } - double sig = 0; - for (T it = beg, ip = p; it != end; ++it, ++ip) { - sig += *ip / ps * (*it - av) * (*it - av); - } - sig = sqrt(sig); - double lower = av - sig * sigs, upper = av + sig * sigs; - double ret = 0, retn = 0; - for (T it = beg, ip = p; it != end; ++it, ++ip) { - if (lower <= *it && *it <= upper) { - ret += *it * *ip; - retn += *ip; - } - } - if (retn <= 1e-10) return av; - return ret / retn; -} - -/*! - * @brief 就只是個取絕對值 - */ -template<class T> -inline T tAbs(T const& t) { - return (t < 0 ? -t : t); -} - -} // meow - -#endif // math_utility_H__ |