add math

42 files changed, 1308 insertions, 321 deletions

diff --git a/meowpp/Usage.h b/meowpp/Usage.h
index 60be936..801b11c 100644
--- a/meowpp/Usage.h
+++ b/meowpp/Usage.h
@@ -1,5 +1,7 @@
-#ifndef   USAGE_H_
-#define   USAGE_H_
+#ifndef   Usage_H__
+#define   Usage_H__
+
+#include <cstdlib>
 
 #include <string>
 #include <vector>
@@ -155,4 +157,4 @@ namespace meow{
 
 #include "Usage.hpp"
 
-#endif // USAGE_H_
+#endif // Usage_H__
diff --git a/meowpp/Usage.hpp b/meowpp/Usage.hpp
index 02f912b..ddf56dc 100644
--- a/meowpp/Usage.hpp
+++ b/meowpp/Usage.hpp
@@ -1,10 +1,14 @@
+#include "Usage.h"
+
+
+#include "utility.h"
+
+#include <cstdlib>
+
+#include <algorithm>
 #include <string>
 #include <vector>
 #include <map>
-#include <algorithm>
-#include <cstdlib>
-
-#include "utility.h"
 
 extern "C"{
 #include <unistd.h>
diff --git a/meowpp/colors/HSL.h b/meowpp/colors/HSL.h
index 2943c1c..db31f0a 100644
--- a/meowpp/colors/HSL.h
+++ b/meowpp/colors/HSL.h
@@ -1,9 +1,11 @@
-#ifndef   HSL_H_
-#define   HSL_H_
+#ifndef   colors_HSL_H_
+#define   colors_HSL_H_
 
 #include "RGB.h"
 #include "YUV.h"
 
+#include <cstdlib>
+
 namespace meow{
   template<class T>
   class HSL{
@@ -63,4 +65,4 @@ namespace meow{
 
 #include "HSL.hpp"
 
-#endif // HSL_H_
+#endif // colors_HSL_H_
diff --git a/meowpp/colors/HSL.hpp b/meowpp/colors/HSL.hpp
index f4b8500..bd9f469 100644
--- a/meowpp/colors/HSL.hpp
+++ b/meowpp/colors/HSL.hpp
@@ -3,7 +3,11 @@
 #include "RGB.h"
 #include "YUV.h"
 
-#include "../utility.h"
+#include "../math/utility.h"
+
+#include <cstdlib>
+
+#include <algorithm>
 
 namespace meow{
   template<class T>
diff --git a/meowpp/colors/HSV.h b/meowpp/colors/HSV.h
index c28cb83..fe3c70d 100644
--- a/meowpp/colors/HSV.h
+++ b/meowpp/colors/HSV.h
@@ -1,10 +1,12 @@
-#ifndef   HSV_H_
-#define   HSV_H_
+#ifndef   colors_HSV_H_
+#define   colors_HSV_H_
 
 #include "RGB.h"
 #include "YUV.h"
 #include "HSL.h"
 
+#include <cstdlib>
+
 namespace meow{
   template<class T>
   class HSV{
@@ -68,4 +70,4 @@ namespace meow{
 
 #include "HSV.hpp"
 
-#endif // HSV_H_
+#endif // colors_HSV_H_
diff --git a/meowpp/colors/HSV.hpp b/meowpp/colors/HSV.hpp
index b028c6c..1838d0d 100644
--- a/meowpp/colors/HSV.hpp
+++ b/meowpp/colors/HSV.hpp
@@ -4,7 +4,11 @@
 #include "YUV.h"
 #include "HSL.h"
 
-#include "../utility.h"
+#include "../math/utility.h"
+
+#include <cstdlib>
+
+#include <algorithm>
 
 namespace meow{
   template<class T>
diff --git a/meowpp/colors/RGB.h b/meowpp/colors/RGB.h
index 5531d96..c7c03c1 100644
--- a/meowpp/colors/RGB.h
+++ b/meowpp/colors/RGB.h
@@ -1,5 +1,7 @@
-#ifndef   RGB_H_
-#define   RGB_H_
+#ifndef   colors_RGB_H_
+#define   colors_RGB_H_
+
+#include <cstdlib>
 
 namespace meow{
   template<class T>
@@ -66,4 +68,4 @@ namespace meow{
 
 #include "RGB.hpp"
 
-#endif // RGB_H_
+#endif // colors_RGB_H_
diff --git a/meowpp/colors/RGB.hpp b/meowpp/colors/RGB.hpp
index f61b7fc..1364814 100644
--- a/meowpp/colors/RGB.hpp
+++ b/meowpp/colors/RGB.hpp
@@ -1,5 +1,11 @@
+#include "RGB.h"
+
+#include <cstdlib>
+
 #include <algorithm>
 
+#include "../math/utility.h"
+
 namespace meow{
   template<class T>
   inline RGB<T>::RGB(){ }
diff --git a/meowpp/colors/YUV.h b/meowpp/colors/YUV.h
index 1b1d60e..0b48375 100644
--- a/meowpp/colors/YUV.h
+++ b/meowpp/colors/YUV.h
@@ -1,8 +1,10 @@
-#ifndef   YUV_H_
-#define   YUV_H_
+#ifndef   colors_YUV_H_
+#define   colors_YUV_H_
 
 #include "RGB.h"
 
+#include <cstdlib>
+
 namespace meow{
   template<class T>
   class YUV{
@@ -58,4 +60,4 @@ namespace meow{
 
 #include "YUV.hpp"
 
-#endif // YUV_H_
+#endif // colors_YUV_H_
diff --git a/meowpp/colors/YUV.hpp b/meowpp/colors/YUV.hpp
index 6847160..9763d2c 100644
--- a/meowpp/colors/YUV.hpp
+++ b/meowpp/colors/YUV.hpp
@@ -1,6 +1,12 @@
-#include <algorithm>
+#include "YUV.h"
+
 #include "RGB.h"
-#include "../utility.h"
+
+#include "../math/utility.h"
+
+#include <cstdlib>
+
+#include <algorithm>
 
 namespace meow{
   template<class T>
diff --git a/meowpp/dsa/BinaryIndexTree.h b/meowpp/dsa/BinaryIndexTree.h
index 4b1b23a..36f24d7 100644
--- a/meowpp/dsa/BinaryIndexTree.h
+++ b/meowpp/dsa/BinaryIndexTree.h
@@ -1,5 +1,9 @@
-#ifndef   __BinaryIndexTree_H__
-#define   __BinaryIndexTree_H__
+#ifndef   dsa_BinaryIndexTree_H__
+#define   dsa_BinaryIndexTree_H__
+
+#include <cstdlib>
+
+#include <vector>
 
 namespace meow{
   //#
@@ -65,5 +69,5 @@ namespace meow{
 
 #include "BinaryIndexTree.hpp"
 
-#endif // BinaryIndexTree_H__
+#endif // dsa_BinaryIndexTree_H__
 
diff --git a/meowpp/dsa/BinaryIndexTree.hpp b/meowpp/dsa/BinaryIndexTree.hpp
index e7e146a..f84a931 100644
--- a/meowpp/dsa/BinaryIndexTree.hpp
+++ b/meowpp/dsa/BinaryIndexTree.hpp
@@ -1,3 +1,10 @@
+#include "BinaryIndexTree.h"
+
+#include <cstdlib>
+
+#include <vector>
+#include <algorithm>
+
 
 namespace meow{
   template<class Value>
diff --git a/meowpp/dsa/DisjointSet.h b/meowpp/dsa/DisjointSet.h
index b918adb..bb777e1 100644
--- a/meowpp/dsa/DisjointSet.h
+++ b/meowpp/dsa/DisjointSet.h
@@ -1,9 +1,9 @@
-#ifndef   DisjointSet_H__
-#define   DisjointSet_H__
+#ifndef   dsa_DisjointSet_H__
+#define   dsa_DisjointSet_H__
 
+#include <cstdlib>
 
 #include <vector>
-#include <cstdlib>
 
 namespace meow{
   //#
@@ -70,4 +70,4 @@ namespace meow{
 #include "DisjointSet.hpp"
 
 
-#endif // DisjointSet_H__
+#endif // dsa_DisjointSet_H__
diff --git a/meowpp/dsa/DisjointSet.hpp b/meowpp/dsa/DisjointSet.hpp
index 3b66626..98b2b98 100644
--- a/meowpp/dsa/DisjointSet.hpp
+++ b/meowpp/dsa/DisjointSet.hpp
@@ -1,3 +1,6 @@
+#include "DisjointSet.h"
+
+
 #include <vector>
 #include <cstdlib>
 
diff --git a/meowpp/dsa/KD_Tree.h b/meowpp/dsa/KD_Tree.h
index 216d9c9..b4e579c 100644
--- a/meowpp/dsa/KD_Tree.h
+++ b/meowpp/dsa/KD_Tree.h
@@ -1,10 +1,10 @@
-#ifndef   KD_Tree_H__
-#define   KD_Tree_H__
+#ifndef   dsa_KD_Tree_H__
+#define   dsa_KD_Tree_H__
 
-#include <vector>
 #include <cstdlib>
+
+#include <vector>
 #include <queue>
-#include "../utility.h"
 
 namespace meow{
   //#
@@ -159,4 +159,4 @@ namespace meow{
 
 #include "KD_Tree.hpp"
 
-#endif // KD_Tree_H__
+#endif // dsa_KD_Tree_H__
diff --git a/meowpp/dsa/KD_Tree.hpp b/meowpp/dsa/KD_Tree.hpp
index 7fea6da..735d4af 100644
--- a/meowpp/dsa/KD_Tree.hpp
+++ b/meowpp/dsa/KD_Tree.hpp
@@ -1,8 +1,14 @@
+#include "KD_Tree.h"
+
+
+#include "../utility.h"
+#include "../math/utility.h"
+
 #include <cstdlib>
+
 #include <vector>
 #include <algorithm>
 #include <queue>
-#include "../utility.h"
 
 namespace meow{
   ////////////////////////////////////////////////////////////////////
diff --git a/meowpp/dsa/MergeableHeap.h b/meowpp/dsa/MergeableHeap.h
index ddcd8a5..0ddb100 100644
--- a/meowpp/dsa/MergeableHeap.h
+++ b/meowpp/dsa/MergeableHeap.h
@@ -1,5 +1,5 @@
-#ifndef   Heap_H__
-#define   Heap_H__
+#ifndef   dsa_MergeableHeap_H__
+#define   dsa_MergeableHeap_H__
 
 #include <cstdlib>
 
@@ -105,4 +105,4 @@ namespace meow{
 
 #include "MergeableHeap.hpp"
 
-#endif // Heap_H__
+#endif // dsa_MergeableHeap_H__
diff --git a/meowpp/dsa/MergeableHeap.hpp b/meowpp/dsa/MergeableHeap.hpp
index e7f5978..1470ac3 100644
--- a/meowpp/dsa/MergeableHeap.hpp
+++ b/meowpp/dsa/MergeableHeap.hpp
@@ -1,5 +1,9 @@
+#include "MergeableHeap.h"
+
 #include <cstdlib>
 
+#include <algorithm>
+
 namespace meow{
   //////////////////////////////////////////////////////////
   //      **# MergeableHeap--Node-- constructor #**       //
diff --git a/meowpp/dsa/SegmentTree.h b/meowpp/dsa/SegmentTree.h
index 4bed51b..52d5a6f 100644
--- a/meowpp/dsa/SegmentTree.h
+++ b/meowpp/dsa/SegmentTree.h
@@ -1,5 +1,9 @@
-#ifndef   SegmentTree_H__
-#define   SegmentTree_H__
+#ifndef   dsa_SegmentTree_H__
+#define   dsa_SegmentTree_H__
+
+#include <cstdlib>
+
+#include <vector>
 
 namespace meow{
   //#
@@ -83,28 +87,8 @@ namespace meow{
       //#|將區間 `[first,last]` 全部都加上 `delta`
       void  offset  (ssize_t __first, ssize_t __last, Value const& __delta);
       
-      //
-      void print(){
-        for(int i = 0; i < _size; i++){
-          query(i, i);
-        }
-        printf("\n");
-        for(int depth = 0, count = 1, size = _size, id = 0;
-            size > 0;
-            depth++, count *= 2, size /= 2){
-          for(int j = 0; j < count; j++, id++){
-            printf("[%3d%c%3d%*s]",
-                   _nodes[id]._value.value,
-                   _nodes[id]._sameFlag ? 's' : ' ',
-                   _nodes[id]._offset.value,
-                   size * 9 - 9,
-                   "");
-          }
-          printf("\n");
-        }
-      }
       
-    //#|=====================================================================
+      //#|=====================================================================
   };
   //#
   //# '''
@@ -113,4 +97,4 @@ namespace meow{
 
 #include "SegmentTree.hpp"
 
-#endif
+#endif // dsa_SegmentTree_H__
diff --git a/meowpp/dsa/SegmentTree.hpp b/meowpp/dsa/SegmentTree.hpp
index a2b74e8..bdd43fb 100644
--- a/meowpp/dsa/SegmentTree.hpp
+++ b/meowpp/dsa/SegmentTree.hpp
@@ -1,5 +1,11 @@
-#include "../utility.h"
+#include "SegmentTree.h"
 
+
+#include "../math/utility.h"
+
+#include <cstdlib>
+
+#include <vector>
 #include <algorithm>
 
 namespace meow{
diff --git a/meowpp/dsa/SplayTree.h b/meowpp/dsa/SplayTree.h
index 38cbe7f..da451b0 100644
--- a/meowpp/dsa/SplayTree.h
+++ b/meowpp/dsa/SplayTree.h
@@ -1,7 +1,9 @@
-#ifndef   SplayTree_h__
-#define   SplayTree_h__
+#ifndef   dsa_SplayTree_h__
+#define   dsa_SplayTree_h__
 
-#include "../utility.h"
+#include <cstdlib>
+
+#include <utility>
 
 namespace meow{
   //#
@@ -53,8 +55,6 @@ namespace meow{
       //
       void  split(Node* __root, Node** __left, Node** __right);
       Node* merge(              Node*  __left, Node*  __right);
-      //
-      void print(Node* __now, int __depth) const;
     public:
       //#==== Custom Type Definitions
       //#
@@ -215,7 +215,6 @@ namespace meow{
       bool merge(SplayTree* __tree2);
       
       
-      void print() const;
       //#|=====================================================================
   };
   //#
@@ -233,4 +232,4 @@ namespace meow{
 
 #include "SplayTree.hpp"
 
-#endif // SplayTree_h__
+#endif // dsa_SplayTree_h__
diff --git a/meowpp/dsa/SplayTree.hpp b/meowpp/dsa/SplayTree.hpp
index de08f63..3b08c14 100644
--- a/meowpp/dsa/SplayTree.hpp
+++ b/meowpp/dsa/SplayTree.hpp
@@ -1,3 +1,9 @@
+#include "SplayTree.h"
+
+
+#include <cstdlib>
+
+#include <utility>
 
 namespace meow{
   ///////////////////////////// **# Node #** /////////////////////////
@@ -161,20 +167,6 @@ namespace meow{
     __left->syncUp();
     return __left;
   }
-  //
-  template<class Key, class Value>
-  inline void SplayTree<Key, Value>::print(Node* __now, int __depth) const{
-    if(__now == NULL) return ;
-    printf("%*s [%llX]:(%lu)\tParent=%llX Left=%llX Right=%llX\n",
-           __depth * 2, "",
-           (long long unsigned)__now,
-           __now->_size,
-           (long long unsigned)__now->_parent,
-           (long long unsigned)__now->_child[0],
-           (long long unsigned)__now->_child[1]);
-    print(__now->_child[0], __depth + 1);
-    print(__now->_child[1], __depth + 1);
-  }
   ///////////////////////// **# Element ##* //////////////////////////
   template<class Key, class Value>
   inline void SplayTree<Key, Value>::Element::reset(Node* __node){
@@ -441,10 +433,5 @@ namespace meow{
     __tree2->_root = NULL;
     return true;
   }
-  template<class Key, class Value>
-  inline void
-  SplayTree<Key, Value>::print() const{
-    print(_root, 0);
-  }
 }
 
diff --git a/meowpp/dsa/SplayTree_Range.h b/meowpp/dsa/SplayTree_Range.h
index e5124c8..0abdc24 100644
--- a/meowpp/dsa/SplayTree_Range.h
+++ b/meowpp/dsa/SplayTree_Range.h
@@ -1,7 +1,11 @@
-#ifndef   SplayTree_Range_h__
-#define   SplayTree_Range_h__
+#ifndef   dsa_SplayTree_Range_h__
+#define   dsa_SplayTree_Range_h__
 
-#include "../utility.h"
+#include <cstdlib>
+
+#include <utility>
+
+#include "../math/utility.h"
 
 namespace meow{
   //#
@@ -57,8 +61,6 @@ namespace meow{
       //
       void  split(Node* __root, Node** __left, Node** __right);
       Node* merge(              Node*  __left, Node*  __right);
-      //
-      void print(Node* __now, int __depth) const;
     public:
       //#==== Custom Type Definitions
       //#
@@ -239,7 +241,6 @@ namespace meow{
       bool merge(SplayTree_Range* __tree2);
 
 
-      void print() const;
       //#|=====================================================================
   };
   //#
@@ -257,4 +258,4 @@ namespace meow{
 
 #include "SplayTree_Range.hpp"
 
-#endif // SplayTree_Range_h__
+#endif // dsa_SplayTree_Range_h__
diff --git a/meowpp/dsa/SplayTree_Range.hpp b/meowpp/dsa/SplayTree_Range.hpp
index 1f216cf..def7ef7 100644
--- a/meowpp/dsa/SplayTree_Range.hpp
+++ b/meowpp/dsa/SplayTree_Range.hpp
@@ -1,3 +1,11 @@
+#include "SplayTree_Range.h"
+
+
+#include <cstdlib>
+
+#include <utility>
+
+#include "../math/utility.h"
 
 namespace meow{
   ///////////////////////////// **# Node #** /////////////////////////
@@ -192,21 +200,6 @@ namespace meow{
     __left->syncUp();
     return __left;
   }
-  //
-  template<class Key, class Value>
-  inline void SplayTree_Range<Key, Value>::print(Node* __now,
-                                                 int __depth) const{
-    if(__now == NULL) return ;
-    printf("%*s [%llX]:(%lu)\tParent=%llX Left=%llX Right=%llX\n",
-           __depth * 2, "",
-           (long long unsigned)__now,
-           __now->_size,
-           (long long unsigned)__now->_parent,
-           (long long unsigned)__now->_child[0],
-           (long long unsigned)__now->_child[1]);
-    print(__now->_child[0], __depth + 1);
-    print(__now->_child[1], __depth + 1);
-  }
   ///////////////////////// **# Element ##* //////////////////////////
   template<class Key, class Value>
   inline void SplayTree_Range<Key, Value>::Element::reset(Node* __node){
@@ -509,10 +502,5 @@ namespace meow{
     __tree2->_root = NULL;
     return true;
   }
-  template<class Key, class Value>
-  inline void
-  SplayTree_Range<Key, Value>::print() const{
-    print(_root, 0);
-  }
 }
 
diff --git a/meowpp/dsa/VP_Tree.h b/meowpp/dsa/VP_Tree.h
index 1877b74..e2daa28 100644
--- a/meowpp/dsa/VP_Tree.h
+++ b/meowpp/dsa/VP_Tree.h
@@ -1,11 +1,13 @@
-#ifndef   VP_Tree_H__
-#define   VP_Tree_H__
+#ifndef   dsa_VP_Tree_H__
+#define   dsa_VP_Tree_H__
+
+#include "../math/utility.h"
+
+#include <cstdlib>
 
 #include <list>
 #include <vector>
-#include <cstdlib>
 #include <queue>
-#include "../utility.h"
 
 namespace meow{
   //#
@@ -90,9 +92,6 @@ namespace meow{
                  Answers*             __out) const;
       void clear(Node* __root);
       Node* dup(Node* __root);
-      //
-      void print(Node* __node, int depth = 1,
-                 Node* __parent = NULL, bool __near = true) const;
     public:
       VP_Tree();
       VP_Tree(VP_Tree const& __tree2);
@@ -152,8 +151,6 @@ namespace meow{
 
 
       //#|=====================================================================
-      
-      void print() const;
   };
   //#
   //#[NOTE]
@@ -170,4 +167,4 @@ namespace meow{
 
 #include "VP_Tree.hpp"
 
-#endif // VP_Tree_H__
+#endif // dsa_VP_Tree_H__
diff --git a/meowpp/dsa/VP_Tree.hpp b/meowpp/dsa/VP_Tree.hpp
index ba97ad7..2026050 100644
--- a/meowpp/dsa/VP_Tree.hpp
+++ b/meowpp/dsa/VP_Tree.hpp
@@ -1,7 +1,10 @@
+#include "VP_Tree.h"
 
 #include <cstdlib>
+
 #include <algorithm>
-#include "../utility.h"
+#include <stack>
+#include "../math/utility.h"
 
 namespace meow{
   ///////////////////// **# Node #** ///////////////////////
@@ -51,7 +54,6 @@ namespace meow{
                                            Scalar const& __b2,
                                            Scalar const& __c2) const{
     // test if  sqrt(__a2) +- sqrt(|__b2|) <= sqrt(__c2)
-    //printf("abc = %lld %lld %lld\n", __a2, __b2, __c2);
     if(__b2 < 0){
       return -distanceCompare(__c2, -__b2, __a2);
     }
@@ -170,29 +172,6 @@ namespace meow{
     ret->_farChild  = dup(__root->_farChild );
     return ret;
   }
-  /////////////////////// **# print #** ////////////////////
-  template<class Vector, class Scalar>
-  inline void
-  VP_Tree<Vector, Scalar>::print(Node* __node, int depth,
-                                 Node* __parent, bool __near) const{
-    if(__node == NULL) return ;
-    printf("%*s%c)Me<%lld,%lld>, rad2 = %7lld",
-           depth * 2, "",
-           __near ? 'N' : 'F',
-           _vectors[__node->_index][0], _vectors[__node->_index][1],
-           __node->_threshold);
-    if(__parent != NULL){
-      printf("   ---<%lld,%lld>:: %lld\n",
-             _vectors[__parent->_index][0],
-             _vectors[__parent->_index][1],
-             distance2(_vectors[__parent->_index],
-                       _vectors[__node  ->_index]));
-    }else{
-      printf("\n");
-    }
-    print(__node->_nearChild, depth + 1, __node, true );
-    print(__node->_farChild , depth + 1, __node, false);
-  }
   ///////// **# construre/destructure/copy oper #** ////////
   template<class Vector, class Scalar>
   inline
@@ -293,13 +272,4 @@ namespace meow{
     _dimension = std::max((size_t)1, __dimension);
     return _dimension;
   }
-  
-  /////////////////////// **# print #** ////////////////////
-  template<class Vector, class Scalar>
-  void inline 
-  VP_Tree<Vector, Scalar>::print() const{
-    printf("\nsize = %lu, dimension = %lu\n", _vectors.size(), _dimension);
-    print(_root, 1);
-    printf("\n\n");
-  }
 };
diff --git a/meowpp/geo/Vector2D.h b/meowpp/geo/Vector2D.h
index 0b18138..7de3523 100644
--- a/meowpp/geo/Vector2D.h
+++ b/meowpp/geo/Vector2D.h
@@ -1,11 +1,10 @@
-#ifndef   Vector2D_H__
-#define   Vector2D_H__
+#ifndef   geo_Vector2D_H__
+#define   geo_Vector2D_H__
 
 #include <cmath>
-#include "../utility.h"
+#include "../math/utility.h"
 
 namespace meow{
-
   template<class Scalar>
   class Vector2D{
     private:
@@ -64,4 +63,4 @@ namespace meow{
 
 }
 
-#endif // Vector2D_H__
+#endif // geo_Vector2D_H__
diff --git a/meowpp/geo/Vector3D.h b/meowpp/geo/Vector3D.h
index f7a669a..08b5705 100644
--- a/meowpp/geo/Vector3D.h
+++ b/meowpp/geo/Vector3D.h
@@ -1,5 +1,5 @@
-#ifndef   Vector3D_H__
-#define   Vector3D_H__
+#ifndef   geo_Vector3D_H__
+#define   geo_Vector3D_H__
 
 #include <cmath>
 #include "../utility.h"
@@ -71,4 +71,4 @@ namespace meow{
   };
 }
 
-#endif // Vector3D_H__
+#endif // geo_Vector3D_H__
diff --git a/meowpp/math/LinearTransformation.h b/meowpp/math/LinearTransformation.h
new file mode 100644
index 0000000..f856f08
--- /dev/null
+++ b/meowpp/math/LinearTransformation.h
@@ -0,0 +1,27 @@
+#ifndef   math_LinearTransformation_H__
+#define   math_LinearTransformation_H__
+
+#include "Transformation.h"
+#include "Matrix.h"
+
+#include <cstdlib>
+
+namespace meow{
+  template<class Scalar>
+  class LinearTransformation: public Transformation<Scalar>{
+    protected:
+      Matrix<Scalar> _matrix;
+      
+      LinearTransformation(size_t  __inputRows, size_t __outputRows,
+                           size_t __psize):
+      Transformation(__inputRows, 1u, __outputRows, 1u, __psize),
+      _matrix(__outputRows, __inputCols, Scalar(0.0)){
+      }
+    public:
+      virtual LinearTransformation(){ }
+      virtual Matrix<Scalar> const&    matrix() const{ return _matrix; }
+      virtual Matrix<Scalar>        invMatrix() const = 0;
+  };
+}
+
+#endif // math_LinearTransformation_H__
diff --git a/meowpp/math/LinearTransformations.h b/meowpp/math/LinearTransformations.h
new file mode 100644
index 0000000..4878834
--- /dev/null
+++ b/meowpp/math/LinearTransformations.h
@@ -0,0 +1,38 @@
+#ifndef   math_LinearTransformations_H__
+#define   math_LinearTransformations_H__
+
+#include "LinearTransformation.h"
+
+#include <cstdlib>
+
+namespace meow{
+  template<class Scalar>
+  class Rotation3D: public LinearTransformation<Scalar>{
+    private:
+      Scalar _theta[3];
+      void calcMatrix();
+    public:
+      Rotation3D();
+      
+      void axisTheta(Matrix<Scalar> const& __axis, Scalar const& __theta);
+      
+      Scalar parameter(size_t __i) const;
+      Scalar parameter(size_t __i, Scalar const& __s);
+      
+      Matrix<Scalar> transformate(Matrix<Scalar> const& __x) const;
+      Matrix<Scalar> jacobian    (Matrix<Scalar> const& __x) const;
+      Matrix<Scalar> jacobian    (Matrix<Scalar> const& __x,
+                                  size_t                __i) const;
+      
+      Matrix<Scalar> invTransformate(Matrix<Scalar> const& __x) const;
+      Matrix<Scalar> invJacobian    (Matrix<Scalar> const& __x) const;
+      Matrix<Scalar> invJacobian    (Matrix<Scalar> const& __x,
+                                     size_t                __i) const;
+      
+      Matrix<Scalar> invMatrix() const;
+  };
+}
+
+#include "LinearTransformations.hpp"
+
+#endif // math_LinearTransformations_H__
diff --git a/meowpp/math/LinearTransformations.hpp b/meowpp/math/LinearTransformations.hpp
new file mode 100644
index 0000000..fcb5b8a
--- /dev/null
+++ b/meowpp/math/LinearTransformations.hpp
@@ -0,0 +1,148 @@
+#include "LinearTransformations.h"
+
+
+#include "utility.h"
+
+#include <cstdlib>
+#include <cmath>
+
+namespace meow{
+  template<class Scalar>
+  inline void
+  Rotation3D<Scalar>::calcMatrix(){
+    Scalar sum(0.0);
+    for(size_t i = 0; i < 3u; i++){
+      sum = sum + squ(_theta(i, 0));
+    }
+    Scalar t(sqrt(double(sum)));
+    std::vector<Scalar> u(3);
+    for(size_t i = 0; i < 3u; i++){
+      u[i] = _theta[i] / t;
+    }
+    Scalar cs(cos(double(t)));
+    Scalar sn(sin(double(t)));
+    
+    _matrix.entry(0, 0, cs + squ(u[0]) * ((1 - cs)));
+    _matrix.entry(1, 1, cs + squ(u[1]) * ((1 - cs)));
+    _matrix.entry(2, 2, cs + squ(u[2]) * ((1 - cs)));
+    _matrix.entry(0, 1, u[0] * u[1] * (1 - cs) - u[2] * sn);
+    _matrix.entry(1, 0, u[1] * u[0] * (1 - cs) + u[2] * sn);
+    _matrix.entry(0, 2, u[0] * u[2] * (1 - cs) + u[1] * sn);
+    _matrix.entry(2, 0, u[2] * u[0] * (1 - cs) - u[1] * sn);
+    _matrix.entry(1, 2, u[1] * u[2] * (1 - cs) - u[0] * sn);
+    _matrix.entry(2, 1, u[2] * u[1] * (1 - cs) + u[0] * sn);
+  }
+  
+  
+  template<class Scalar>
+  inline
+  Rotation3D<Scalar>::Rotation3D():
+  LinearTransformation(3u, 3u, 3u){
+    _theta[0] = _theta[1] = _theta[2] = Scalar(0);
+    _matrix.identitied();
+  }
+  
+  
+  template<class Scalar>
+  inline void
+  Rotation3D<Scalar>::axisTheta(Matrix<Scalar> const& __axis,
+                                Scalar         const& __theta){
+    Scalar sum(0.0);
+    for(size_t i = 0; i < 3u; i++)
+      sum = sum + squ(__axis(i, 0));
+    Scalar t(sqrt(double(sum)));
+    for(size_t i = 0; i < 3u; i++)
+      _theta[i] = __axis(i, 0) * __theta / t;
+    calcMatrix();
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  Rotation3D<Scalar>::parameter(size_t __i) const{
+    return _theta[__i];
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  Rotation3D<Scalar>::parameter(size_t __i, Scalar const& __s) const{
+    _theta[__i] =  __s;
+    calcMatrix();
+    return _theta[__i];
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  Rotation3D<Scalar>::transformate(Matrix<Scalar> const& __x) const{
+    return _matrix * __x;
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  Rotation3D<Scalar>::jacobian(Matrix<Scalar> const& __x) const{
+    return _matrix;
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  Rotation3D<Scalar>::jacobian(Matrix<Scalar> const& __x,
+                               size_t                __i) const{
+    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 _matrix * mid * __x;
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  Rotation3D<Scalar>::invTransformate(Matrix<Scalar> const& __x) const{
+    return _matrix.transpose() * __x;
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  Rotation3D<Scalar>::invJacobian(Matrix<Scalar> const& __x) const{
+    return _matrix.transpose();
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  Rotation3D<Scalar>::invJacobian(Matrix<Scalar> const& __x,
+                                  size_t                __i) const{
+    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 _matrix.transpose() * mid.transpose() * __x;
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  Rotation3D<Scalar>::invMatrix() const{
+    return _matrix.transpose();
+  }
+}
+
diff --git a/meowpp/math/Matrix.h b/meowpp/math/Matrix.h
new file mode 100644
index 0000000..0a52d8c
--- /dev/null
+++ b/meowpp/math/Matrix.h
@@ -0,0 +1,74 @@
+#ifndef   math_Matrix_H__
+#define   math_Matrix_H__
+
+
+#include <cstdlib>
+
+#include <vector>
+
+namespace meow{
+  template<class Entry>
+  class Matrix{
+    private:
+      size_t             _rows;
+      size_t             _cols;
+      std::vector<Entry> _entries;
+      //
+      size_t index(size_t __r, size_t __c) const;
+    public:
+      Matrix();
+      Matrix(Matrix const& __m);
+      Matrix(size_t __rows, size_t __cols, Entry const& __entry);
+      ~Matrix();
+      
+      Matrix& copy(Matrix const& __m);
+      
+      bool valid() const;
+      
+      size_t rows() const;
+      size_t cols() const;
+      size_t size(size_t __rows, size_t __cols);
+      
+      Entry const& entry(size_t __i, size_t __j) const;
+      Entry const& entry(size_t __i, size_t __j, Entry const& __entry);
+      void entries(size_t __rFirst, size_t __rLast,
+                   size_t __cFirst, size_t __cLast,
+                   Entry const& __entry);
+      
+      Matrix subMatrix(size_t __rFirst, size_t __rLast,
+                       size_t __cFirst, size_t __cLast) const;
+      Matrix row(size_t __row) const;
+      Matrix col(size_t __col) const;
+      
+      Matrix positive() const;
+      Matrix negative() const;
+      
+      Matrix add(Matrix const& __m) const;
+      Matrix sub(Matrix const& __m) const;
+      Matrix mul(Matrix const& __m) const;
+      
+      Matrix mul(Entry const& __s) const;
+      Matrix div(Entry const& __s) const;
+      
+      Matrix  identity  () const;
+      Matrix& identitied();
+      
+      Matrix  transpose () const;
+      Matrix& transposed();
+      
+      Matrix& operator=(Matrix const& __m);
+      Entry const& operator()(size_t __i, size_t __j) const;
+      Entry      & operator()(size_t __i, size_t __j);
+      Matrix operator+() const;
+      Matrix operator-() const;
+      Matrix operator+(Matrix const& __m) const;
+      Matrix operator-(Matrix const& __m) const;
+      Matrix operator*(Matrix const& __m) const;
+      Matrix operator*(Entry const& __s) const;
+      Matrix operator/(Entry const& __s) const;
+  };
+}
+
+#include "Matrix.hpp"
+
+#endif // math_Matrix_H__
diff --git a/meowpp/math/Matrix.hpp b/meowpp/math/Matrix.hpp
new file mode 100644
index 0000000..8de3057
--- /dev/null
+++ b/meowpp/math/Matrix.hpp
@@ -0,0 +1,329 @@
+#include "Matrix.h"
+
+#include <cstdlib>
+
+namespace meow{
+  template<class Entry>
+  inline size_t
+  Matrix<Entry>::index(size_t __r, size_t __c) const{
+    return __r * _cols + __c;
+  }
+  
+  
+  template<class Entry>
+  inline
+  Matrix<Entry>::Matrix():
+  _rows(0),
+  _cols(0),
+  _entries(0){
+  }
+  
+  
+  template<class Entry>
+  inline
+  Matrix<Entry>::Matrix(Matrix const& __m):
+  _rows(__m._rows),
+  _cols(__m._cols),
+  _entries(__m._entries){
+  }
+  
+  
+  template<class Entry>
+  inline
+  Matrix<Entry>::Matrix(size_t __rows, size_t __cols, Entry const& __entry):
+  _rows(__rows),
+  _cols(__cols),
+  _entries(__rows * __cols, __entry){
+  }
+  
+  
+  template<class Entry>
+  inline
+  Matrix<Entry>::~Matrix(){
+  }
+  
+  template<class Entry>
+  inline Matrix<Entry>&
+  Matrix<Entry>::copy(Matrix const& __m){
+    _rows    = __m._rows;
+    _cols    = __m._cols;
+    _entries = __m._entries;
+  }
+  
+  
+  template<class Entry>
+  inline bool
+  Matrix<Entry>::valid() const{
+    return (rows() > 0 && cols() > 0);
+  }
+  
+  template<class Entry>
+  inline size_t
+  Matrix<Entry>::rows() const{
+    return _rows;
+  }
+  template<class Entry>
+  inline size_t
+  Matrix<Entry>::cols() const{
+    return _cols;
+  }
+  
+  
+  template<class Entry>
+  inline size_t
+  Matrix<Entry>::size(size_t __rows, size_t __cols){
+    if(__rows != rows() || __cols != cols()){
+      if(__cols == cols()){
+        _entries.resize(_rows * (_cols = __cols));
+      }else{
+        Matrix new_matrix(*this);
+        size_t R = std::min(rows(), __rows);
+        size_t C = std::min(cols(), __cols);
+        for(size_t r = 0; r < R; r++)
+          for(size_t c = 0; c < C; c++)
+            new_matrix.entry(r, c, entry(r, c));
+        copy(new_matrix);
+      }
+    }
+    return rows() * cols();
+  }
+  
+  
+  template<class Entry>
+  inline Entry const&
+  Matrix<Entry>::entry(size_t __i, size_t __j) const{
+    return _entries[index(__i, __j)];
+  }
+  
+  
+  template<class Entry>
+  inline Entry const&
+  Matrix<Entry>::entry(size_t __i, size_t __j, Entry const& __entry){
+    _entries[index(__i, __j)] = __entry;
+    return entry(__i, __j);
+  }
+
+
+  template<class Entry>
+  inline void
+  Matrix<Entry>::entries(size_t __rFirst, size_t __rLast,
+                         size_t __cFirst, size_t __cLast,
+                         Entry const& __entry){
+    for(size_t r = __rFirst; r <= __rLast; r++)
+      for(size_t c = __cFirst; c <= __cLast; c++)
+        entry(r, c, __entry);
+  }
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::subMatrix(size_t __rFirst, size_t __rLast,
+                           size_t __cFirst, size_t __cLast) const{
+    if(__rFirst > __rLast || __cFirst > __cLast)
+      return Matrix<Entry>();
+    Matrix ret(__rLast - __rFirst + 1, __cLast - __cFirst + 1);
+    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;
+  }
+  
+  
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::row(size_t __row) const{
+    return subMatrix(__row, __row, 0, cols() - 1);
+  }
+  
+  
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::col(size_t __col) const{
+    return subMatrix(0, rows() - 1, __col, __col);
+  }
+
+
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::positive() const{
+    return *this;
+  }
+
+
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::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;
+  }
+
+
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::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;
+  }
+  
+
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::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;
+  }
+  
+  
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::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;
+  }
+
+  
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::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;
+  }
+  
+  
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::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;
+  }
+  
+  
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::identity() const{
+    Matrix ret(rows(), cols(), Entry(0));
+    for(size_t rc = 0, RC = std::min(rows(), cols()); rc < RC; rc++)
+      ret.entry(rc, rc, Entry(1));
+    return ret;
+  }
+  
+  
+  template<class Entry>
+  inline Matrix<Entry>&
+  Matrix<Entry>::identitied(){
+    copy(identity());
+    return (*this);
+  }
+
+  
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::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;
+  }
+  
+  
+  template<class Entry>
+  inline Matrix<Entry>&
+  Matrix<Entry>::transposed(){
+    copy(transpose());
+    return (*this);
+  }
+  
+  
+  template<class Entry>
+  inline Matrix<Entry>&
+  Matrix<Entry>::operator=(Matrix const& __m){
+    return copy(__m);
+  }
+
+  
+  template<class Entry>
+  inline Entry const&
+  Matrix<Entry>::operator()(size_t __i, size_t __j) const{
+    return entry(__i, __j);
+  }
+  
+  
+  template<class Entry>
+  inline Entry&
+  Matrix<Entry>::operator()(size_t __i, size_t __j){
+    return _entries[index(__i, __j)];
+  }
+  
+  
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::operator+() const{
+    return positive();
+  }
+
+
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::operator-() const{
+    return negative();
+  }
+  
+  
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::operator+(Matrix const& __m) const{
+    return add(__m);
+  }
+  
+  
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::operator-(Matrix const& __m) const{
+    return sub(__m);
+  }
+  
+  
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::operator*(Matrix const& __m) const{
+    return mul(__m);
+  }
+  
+  
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::operator*(Entry const& __s) const{
+    return mul(__s);
+  }
+  
+  
+  template<class Entry>
+  inline Matrix<Entry>
+  Matrix<Entry>::operator/(Entry const& __s) const{
+    return div(__s);
+  }
+}
diff --git a/meowpp/math/Transformation.h b/meowpp/math/Transformation.h
new file mode 100644
index 0000000..abf2649
--- /dev/null
+++ b/meowpp/math/Transformation.h
@@ -0,0 +1,59 @@
+#ifndef   math_Transformation_H__
+#define   math_Transformation_H__
+
+#include "Matrix.h"
+
+#include <cstdlib>
+
+namespace meow{
+  template<class Scalar>
+  class Transformation{
+    private:
+      size_t  _inputRows;
+      size_t  _inputCols;
+      size_t _outputRows;
+      size_t _outputCols;
+      size_t _psize;
+    protected:
+      Transformation(size_t  __inputRows, size_t  __inputCols,
+                     size_t __outputRows, size_t __outputCols,
+                     size_t __psize):
+      _inputRows (__inputRows),
+      _inputCols (__inputCols),
+      _outputRows(__outputRows),
+      _outputCols(__outputCols),
+      _psize(__psize){
+      }
+    public:
+      virtual Transformation(){ }
+      
+      size_t     inputRows() const{ return  _inputRows; }
+      size_t     inputCols() const{ return  _inputCols; }
+      size_t    outputRows() const{ return _outputRows; }
+      size_t    outputCols() const{ return _outputCols; }
+      size_t parameterSize() const{ return _psize; }
+      
+      virtual Scalar parameter(size_t __i) const = 0;
+      virtual Scalar parameter(size_t __i, Scalar const& __s) = 0;
+      
+      virtual Matrix<Scalar> transformate(Matrix<Scalar> const& __x) const = 0;
+      virtual Matrix<Scalar> jacobian    (Matrix<Scalar> const& __x) const = 0;
+      virtual Matrix<Scalar> jacobian    (Matrix<Scalar> const& __x,
+                                          size_t                __i) const = 0;
+      
+      virtual bool inversable() const{ return false; }
+      
+      virtual Matrix<Scalar> invTransformate(Matrix<Scalar> const& __x) const{
+        return Matrix<Scalar>();
+      }
+      virtual Matrix<Scalar> invJacobian(Matrix<Scalar> const& __x) const{
+        return Matrix<Scalar>();
+      }
+      virtual Matrix<Scalar> invJacobian(Matrix<Scalar> const& __x,
+                                         size_t                __i) const{
+        return Matrix<Scalar>();
+      }
+  };
+}
+
+#endif // math_Transformation_H__
diff --git a/meowpp/math/Transformations.h b/meowpp/math/Transformations.h
new file mode 100644
index 0000000..3b20ebe
--- /dev/null
+++ b/meowpp/math/Transformations.h
@@ -0,0 +1,53 @@
+#ifndef   math_Transformations_H__
+#define   math_Transformations_H__
+
+#include "Transformation.h"
+
+#include <cstdlib>
+
+namespace meow{
+  template<class Scalar>
+  class BallProjection: public Transformation<Scalar>{
+    private:
+      Scalar _radius;
+      size_t _dimension;
+    public:
+      BallProjection(size_t __dimension);
+      BallProjection(size_t __dimension, Scalar const& __radius);
+      
+      Scalar parameter(size_t __i) const;
+      Scalar parameter(size_t __i, Scalar const& __s);
+      
+      Scalar radius() const;
+      Scalar radius(Scalar const& __radius);
+      
+      Matrix<Scalar> transformate(Matrix<Scalar> const& __x) const;
+      Matrix<Scalar> jacobian    (Matrix<Scalar> const& __x            ) const;
+      Matrix<Scalar> jacobian    (Matrix<Scalar> const& __x, size_t __i) const;
+  };
+  
+  
+  template<class Scalar>
+  class PhotoProjection: public Transformation<Scalar>{
+    private:
+      Scalar _focal;
+      size_t _dimension;
+    public:
+      PhotoProjection(size_t __dimension);
+      PhotoProjection(size_t __dimension, Scalar const& __focal);
+      
+      Scalar parameter(size_t __i) const;
+      Scalar parameter(size_t __i, Scalar const& __s);
+      
+      Scalar focal() const;
+      Scalar focal(Scalar const& __focal);
+      
+      Matrix<Scalar> transformate(Matrix<Scalar> const& __x) const;
+      Matrix<Scalar> jacobian    (Matrix<Scalar> const& __x            ) const;
+      Matrix<Scalar> jacobian    (Matrix<Scalar> const& __x, size_t __i) const;
+  };
+}
+
+#include "Transformations.hpp"
+
+#endif // Transformations_H__
diff --git a/meowpp/math/Transformations.hpp b/meowpp/math/Transformations.hpp
new file mode 100644
index 0000000..645b1a8
--- /dev/null
+++ b/meowpp/math/Transformations.hpp
@@ -0,0 +1,193 @@
+#include "Transformations.h"
+#include "utility.h"
+
+#include <cstdlib>
+
+namespace meow{
+  template<class Scalar>
+  inline
+  BallProjection<Scalar>::BallProjection(size_t __dimension):
+  Transformation(__dimension, 1, __dimension, 1, 1),
+  _dimension(__dimension),
+  _radius(1){
+  }
+  
+  
+  template<class Scalar>
+  inline
+  BallProjection<Scalar>::BallProjection(size_t __dimension,
+                                         Scalar const& __radius):
+  Transformation(__dimension, 1, __dimension, 1, 1),
+  _dimension(__dimension),
+  _radius(1){
+    radius(__radius);
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  BallProjection<Scalar>::parameter(size_t __i) const{
+    return radius();
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  BallProjection<Scalar>::parameter(size_t __i, Scalar const& __s){
+    return radius(__s);
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  BallProjection<Scalar>::radius() const{
+    return _radius;
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  BallProjection<Scalar>::radius(Scalar const& __radius){
+    _radius = __radius;
+    return radius();
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  BallProjection<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 < _dimension; i++){
+        sum = sum + squ(ret(i, c));
+      }
+      Scalar len(sqrt(double(sum)));
+      for(size_t i = 0; i < _dimension; i++){
+        ret(i, c) = ret(i, c) * radius() / len;
+      }
+    }
+    return ret;
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  BallProjection<Scalar>::jacobian(Matrix<Scalar> const& __x) const{
+    Scalar sum(0);
+    for(size_t i = 0; i < _dimension; 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++)
+      for(size_t j = 0; j < _dimension; j++)
+        ret(i, j) = ((i == j)
+                     ? 1.0 / len - squ(__x(i, 0)) / cub(len) * radius()
+                     : - __x(i, 0) * __x(j, 0) / cub(len) * radius()
+                     );
+    return ret;
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  BallProjection<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++){
+      sum = sum + squ(__x(i, 0));
+    }
+    Scalar len(sqrt(double(sum)));
+    for(size_t i = 0; i < _dimension; i++){
+      ret(i, 0) = ret(i, 0) / len;
+    }
+    return ret;
+  }
+  
+  
+  template<class Scalar>
+  inline
+  PhotoProjection<Scalar>::PhotoProjection(size_t __dimension):
+  Transformation(__dimension, 1, __dimension, 1, 0),
+  _dimension(__dimension),
+  _focal(1){
+  }
+  
+  
+  template<class Scalar>
+  inline
+  PhotoProjection<Scalar>::PhotoProjection(size_t        __dimension,
+                                           Scalar const& __focal):
+  Transformation(__dimension, 1, __dimension, 1, 0),
+  _dimension(__dimension),
+  _focal(1){
+    focal(__focal);
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  PhotoProjection<Scalar>::parameter(size_t __i) const{
+    return focal();
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  PhotoProjection<Scalar>::parameter(size_t __i, Scalar const& __s){
+    return focal(__s);
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  PhotoProjection<Scalar>::focal() const{
+    return _focal;
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  PhotoProjection<Scalar>::focal(Scalar const& __focal){
+    _focal = __focal;
+    return focal();
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  PhotoProjection<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++){
+        ret(i, c) = ret(i, c) * _focal / ret(_dimension - 1, c);
+      }
+    }
+    return ret;
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  PhotoProjection<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(_dimension - 1, 0) ;
+      ret(i, _dimension - 1) = -_focal / squ(__x(_dimension - 1, 0));
+    }
+    return ret;
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  PhotoProjection<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++){
+      ret(i, 0) = ret(i, 0) / ret(_dimension - 1, 0);
+    }
+    return ret;
+  }
+}
diff --git a/meowpp/math/utility.h b/meowpp/math/utility.h
new file mode 100644
index 0000000..9bfbb7a
--- /dev/null
+++ b/meowpp/math/utility.h
@@ -0,0 +1,82 @@
+#ifndef   math_utility_H__
+#define   math_utility_H__
+
+#include <cstdlib>
+
+namespace meow{
+  static const double PI = 3.14159265358979323846264338327950288;
+  //#
+  //# === meow:: *Functios* in math/utility.h
+  //#
+  //# [options="header",width="100%",cols="1>s,5<,1<,10<",grid="rows"]
+  //# |==============================================================
+  //# |Name             | Parameters                | Return_Type | Description
+
+  
+  //# |noEPS<T>       |(T `value`, T `eps` = 1e-9) | T |
+  //# 如果abs(輸入的數值) < eps, 則回傳0, 否則回傳輸入的數值
+  template<class T>
+  inline T noEPS(T value, T eps = 1e-9);
+  
+  
+  //# |normalize<T>   |(T `lower`, T `upper`, \ T value)
+  //# | T | `(value - lower) / (upper - lower)`
+  template<class T>
+  inline T normalize(T lower, T upper, T value);
+  
+  
+  //# |denormalize<T> |(T `lower`, T `upper`,
+  //# \ T `ratio`) | T | `lower + (upper - lower) * ratio`
+  template<class T>
+  inline T denormalize(T lower, T upper, T ratio);
+  
+  
+  //# |ratioMapping<T>|(T `l1`, T `u1`,
+  //# \T `m1`, T `l2`,\T `u2`)
+  //# | T | `denormalize(l2, u2, normalize(l1, u1, m1))`
+  template<class T>
+  inline T ratioMapping(T l1, T u1, T m1, T l2, T u2);
+  
+  
+  //# |inRange<T>     |(T const& `mn`, T const& `mx`, \ T const& `v`) | T | 
+  //# `std::max(mn, std::min(mx, v))`
+  template<class T>
+  inline T inRange(T const& mn, T const& mx, T const& v);
+  
+  
+  //# |squ<T>  |(T const& `x`) | T| `x * x`
+  template<class T>
+  inline T squ(T const& x);
+  
+  
+  //# |cub<T>  |(T const& `x`) | T| `x * x * x`
+  template<class T>
+  inline T cub(T const& x);
+
+  //# |average<T>|(T const& `beg`, T const& `end`, \ double `sigs`)| T|
+  //# 只將 `sigs` 個標準差以內的數據拿來取平均
+  template<class T>
+  inline double average(T const& beg, T const& end, double sigs);
+  
+  
+  //# |average<T>|(T const& `beg`, T const& `end`,
+  //# \ T const& `p`, double `sigs`)| T| 同上, 不過這次用 `p` 來加權平均
+  template<class T>
+  inline double average(T const& beg, T const& end, T const& p, double sigs);
+  
+  //# |==============================================================
+  
+  
+  //# 
+  //# [NOTE]
+  //# ====================================
+  //#  * 額外附贈一個 `const double PI = 3.141592653589......`
+  //# ====================================
+  //#
+  //# '''
+  //#
+}
+
+#include "utility.hpp"
+
+#endif // math_utility_H__
diff --git a/meowpp/math/utility.hpp b/meowpp/math/utility.hpp
new file mode 100644
index 0000000..90cbdce
--- /dev/null
+++ b/meowpp/math/utility.hpp
@@ -0,0 +1,102 @@
+#include "utility.h"
+
+#include <algorithm>
+
+#include <cmath>
+
+namespace meow{
+
+  template<class T>
+  inline T noEPS(T value, T eps){
+    T epsp((eps < T(0)) ? -eps : eps);
+    return ((value < -epsp || value > epsp)
+            ? value
+            : T(0));
+  }
+
+  template<class T>
+  inline T normalize(T lower, T upper, T value){
+    return (value - lower) / (upper - lower);
+  }
+
+  template<class T>
+  inline T denormalize(T lower, T upper, T _ratio){
+    return lower + _ratio * (upper - lower);
+  }
+
+  template<class T>
+  inline T ratioMapping(T l1, T u1, T m1, T l2, T u2){
+    return denormalize(l2, u2, normalize(l1, u1, m1));
+  }
+
+  template<class T>
+  inline T inRange(T const& mn, T const& mx, T const& v){
+    return std::min(mx, std::max(mn, v));
+  }
+  
+  
+  template<class T>
+  inline T squ(T const& x){
+    return x * x;
+  }
+  
+  
+  template<class T>
+  inline T cub(T const& x){
+    return x * x * x;
+  }
+  
+  
+  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;
+  }
+  
+  
+  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;
+  }
+}
diff --git a/meowpp/oo/Register_Implement.h b/meowpp/oo/Register_Implement.h
index a27591d..910ee11 100644
--- a/meowpp/oo/Register_Implement.h
+++ b/meowpp/oo/Register_Implement.h
@@ -1,5 +1,5 @@
-#ifndef   REGISTER_IMPLEMENT_H_
-#define   REGISTER_IMPLEMENT_H_
+#ifndef   oo_Register_Implement_H__
+#define   oo_Register_Implement_H__
 
 #include <map>
 #include <vector>
@@ -52,4 +52,4 @@ namespace meow{
 
 #include "Register_Implement.hpp"
 
-#endif // REGISTER_IMPLEMENT_H_
+#endif // oo_Register_Implement_H__
diff --git a/meowpp/oo/Register_Implement.hpp b/meowpp/oo/Register_Implement.hpp
index 206b0a3..bee7f3e 100644
--- a/meowpp/oo/Register_Implement.hpp
+++ b/meowpp/oo/Register_Implement.hpp
@@ -1,3 +1,6 @@
+#include "Register_Implement.h"
+
+
 #include <map>
 #include <vector>
 
diff --git a/meowpp/utility.h b/meowpp/utility.h
index 12909e9..14b7af1 100644
--- a/meowpp/utility.h
+++ b/meowpp/utility.h
@@ -1,110 +1,69 @@
-#ifndef   UTILITY_H_
-#define   UTILITY_H_
+#ifndef   utility_H__
+#define   utility_H__
+
+#include <cstdlib>
 
 #include <string>
-#include <stack>
-#include <cctype>
 
 namespace meow{
-  inline std::string stringPrintf(char const * fmt, ...);
-  inline std::string stringReplace(std::string str,
-                                   std::string const& from,
-                                   std::string const& to);
-  inline bool cstringEndWith(char const* str, int n, ...);
-#define debugPrintf(...) \
-  meow::debugPrintf_(\
-                     __FILE__,\
-                     __PRETTY_FUNCTION__,\
-                     __LINE__,\
-                     meow::stringPrintf(__VA_ARGS__).c_str())
-  inline void debugPrintf_(char const* file,
-                           char const* func,
-                           size_t      line,
-                           char const* msg);
-  inline void messagePrintf(int level_change, char const* fmt, ...);
-
-  static const double PI = 3.14159265358979323846264338327950288;
-
-  inline double noEPS(double value, double eps = 1e-9);
-  inline double   normalize(double lower, double upper, double value);
-  inline double denormalize(double lower, double upper, double ratio);
-  inline double ratioMapping(double l1, double u1, double m1,
-                             double l2, double u2);
-  template<class T>
-  inline T inRange(T const& mn, T const& mx, T const& v);
-  template<class T>
-  inline T squ(T const& x);
-
-  template<class T>
-  inline double average(T const& beg, T const& end, double sigs);
-  template<class T>
-  inline double average(T const& beg, T const& end, T const& p, double sigs);
-
   //#
   //# === meow:: *Functios* in utility.h
   //#
   //# [options="header",width="100%",cols="1>s,5<,1<,10<",grid="rows"]
   //# |==============================================================
   //# |Name             | Parameters                | Return_Type | Description
+  
+  
   //# |stringPrintf     |(char const * `fmt`, ...)  | std::string
   //# |Format print to C++ string and return it
+  inline std::string stringPrintf(char const * fmt, ...);
   
   
   //# |stringReplace  |(std::string `str`,\std::string const&
   //#                   `from`,\std::string const& `to`)    | std::string
   //# |Return a string like `str`, but all `from` be replaced by `to`
+  inline std::string stringReplace(std::string str,
+                                   std::string const& from,
+                                   std::string const& to);
   
   
   //# |cstringEndWith |(char const* `str`,\int `n`, ...) | bool 
   //# |Return whether `str` is end with one of the c-string you specify in
   //#  the parameters or not
+  inline bool cstringEndWith(char const* str, int n, ...);
   
   
   //# |debugPrintf    |(char const* `fmt`, ...)     | void
   //# |Print debug message (file name, line number, ...etc) when `DEBUG` is
   //#  defined
+#define debugPrintf(...) \
+  meow::debugPrintf_(\
+                     __FILE__,\
+                     __PRETTY_FUNCTION__,\
+                     __LINE__,\
+                     meow::stringPrintf(__VA_ARGS__).c_str())
+  inline void debugPrintf_(char const* file,
+                           char const* func,
+                           size_t      line,
+                           char const* msg);
   
   
   //# |messagePrintf  |(int `level_change`,\char const* `fmt`, ...) | void
   //# |階層式的訊息輸出
+  inline void messagePrintf(int level_change, char const* fmt, ...);
+
   
-  
-  //# |noEPS          |(double `value`, double `eps` = 1e-9) | double |
-  //# 如果abs(輸入的數值) < eps, 則回傳0, 否則回傳輸入的數值
-  
-  
-  //# |normalize      |(double `lower`, double `upper`, \ double value)
-  //# | double | `(value - lower) / (upper - lower)`
-  
-  
-  //# |denormalize    |(double `lower`, double `upper`,
-  //# \ double `ratio`) | double | `lower + (upper - lower) * ratio`
-  
-  
-  //# |ratioMapping   |(double `l1`, double `u1`,
-  //# \double `m1`, double `l2`,\double `u2`)
-  //# | double | `denormalize(l2, u2, normalize(l1, u1, m1))`
-  
-  
-  //# |inRange<T>     |(T const& `mn`, T const& `mx`, \ T const& `v`) | T | 
-  //# `std::max(mn, std::min(mx, v))`
-  
-  
-  //# |squ<T>  |(T const& `x`) | T| `x * x`
-  
-  
-  //# |average<T>|(T const& `beg`, T const& `end`, \ double `sigs`)| T|
-  //# 只將 `sigs` 個標準差以內的數據拿來取平均
+  //# |filenameCompare |(std::string const& `f1`, std::string const& `f2`)|void
+  //# | 依照 `a0.txt < a1.txt < a2.txt < a10.txt` 的字串比較方法比較字串
+  inline bool filenameCompare(std::string const& f1, std::string const& f2);
   
   
-  //# |average<T>|(T const& `beg`, T const& `end`,
-  //# \ T const& `p`, double `sigs`)| T| 同上, 不過這次用 `p` 來加權平均
   //# |==============================================================
+  
   //# 
   //# [NOTE]
   //# ====================================
   //#  * `stringReplace()` 不是用什麼好方法寫的因此執行效率很低請別虐待它. 
-  //#  * 額外附贈一個 `const double PI = 3.141592653589......`
   //# ====================================
   //#
   //# '''
@@ -113,4 +72,4 @@ namespace meow{
 
 #include "utility.hpp"
 
-#endif // UTILITY_H_
+#endif // utility_H__
diff --git a/meowpp/utility.hpp b/meowpp/utility.hpp
index 9f6d708..141de41 100644
--- a/meowpp/utility.hpp
+++ b/meowpp/utility.hpp
@@ -1,11 +1,13 @@
-#include <string>
-#include <stack>
+#include "utility.h"
+
+
 #include <cstdio>
+#include <cstdlib>
 #include <cstdarg>
-#include <algorithm>
-#include <cctype>
 #include <cstring>
-#include <cmath>
+
+#include <string>
+#include <stack>
 
 namespace meow{
 
@@ -18,6 +20,7 @@ namespace meow{
     return std::string(str);
   }
 
+
   inline std::string stringReplace(std::string str,
                                    std::string const& from,
                                    std::string const& to){
@@ -29,6 +32,7 @@ namespace meow{
     return out;
   }
 
+
   inline bool cstringEndWith(char const* str, int n, ...){
     int len = strlen(str);
     va_list args;
@@ -44,6 +48,7 @@ namespace meow{
     return false;
   }
 
+
   inline void debugPrintf_(char const* file,
                            char const* func,
                            size_t      line,
@@ -53,6 +58,7 @@ namespace meow{
 #endif // DEBUG
   }
 
+
   inline void messagePrintf(int level_change, char const* fmt, ...){
     static int      level =  0;
     static int last_level = -5;
@@ -75,23 +81,6 @@ namespace meow{
     fflush(stdout);
   }
 
-  inline double noEPS(double value, double eps){
-    return (fabs(value) <= fabs(eps) ? 0 : value);
-  }
-
-  inline double normalize(double lower, double upper, double value){
-    return (value - lower) / (upper - lower);
-  }
-
-  inline double denormalize(double lower, double upper, double ratio){
-    return lower + ratio * (upper - lower);
-  }
-
-  inline double ratioMapping(double l1, double u1, double m1,
-                             double l2, double u2){
-    return denormalize(l2, u2, normalize(l1, u1, m1));
-  }
-
   inline bool filenameCompare(std::string const& f1, std::string const& f2){
     char const* s1 = f1.c_str();
     char const* s2 = f2.c_str();
@@ -115,62 +104,4 @@ namespace meow{
     }
     return false;
   }
-  template<class T>
-  inline T inRange(T const& mn, T const& mx, T const& v){
-    return std::min(mx, std::max(mn, v));
-  }
-  template<class T>
-  inline T squ(T const& x){
-    return x * x;
-  }
-  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;
-  }
-  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;
-  }
 }