x

41 files changed, 1779 insertions, 1148 deletions

diff --git a/meowpp/Self.h b/meowpp/Self.h
index 5a62984..7de568b 100644
--- a/meowpp/Self.h
+++ b/meowpp/Self.h
@@ -15,32 +15,23 @@ namespace meow {
  * private:
  *   struct Myself {
  *     int data;
- *     Myself(){
+ *     Myself() {                               // 必要
  *       data = 0;
  *     }
+ *     Myself(Myself const& b): data(b.data) {  // 必要, copy constructor
+ *     }
  *     ~Myself() {
  *     }
- *     Myself copyFrom(Myself const& ms) const {
- *       data = ms.data;
+ *     bool operator==(Myself const& b) const { // 不一定需要(有用到equal()才要)
+ *       return (data == b.data);
  *     }
  *   };
+ *
  *   Self<Myself> const self;
  * public:
- *   A(): self(true) { // self(true) 表示要建立實體, 即struct Myself
- *   }
- *   A(A const& a): self(false) { // for std::swap
- *     copyFrom(a);
- *   }
- *   // A(A const& a); // disable 模糊用法, 與上者二選一
- *   A(A const& a, bool reference): self(false) {
- *     if (reference) {
- *       referenceFrom(a);
- *     } else {
- *       copyFrom(a);
- *     }
- *   }
- *   ~A() {
- *   }
+ *   A(): self() { }                                // 預設呼叫Myself()
+ *   A(A const& a): self(a.self, COPY_FROM) { }     // 用Myself(b)做到copyFrom
+ *   ~A() { }
  *   void setMemeber(int k) {
  *     self()->data = k;  // self()->?? 可以有write權限
  *   }
@@ -56,138 +47,189 @@ namespace meow {
  *   A& operator=(A const& b) { // for std::swap
  *     copyFrom(b);
  *   }
- *   A& operator=(A const& b); // 避免諢亂用法
  * };
  * @endcode
  *
- *@author cathook
+ * @author cathook
  *
- *@warning \c Self 這個class會把\c operator= 給disable掉, 所以使用它當
- *         kernel的class預設的 \c operator= 都會無法使用
+ * @warning \c Self 這個class會把\c operator= 給disable掉, 所以使用它當
+ *          kernel的class預設的 \c operator= 都會無法使用
  */
 template<class Data>
 class Self {
+public:
+  /*!
+   * @brief 複製資料的方法
+   */
+  enum DuplicateType {
+    COPY_FROM,          //!< 跟一般的複製一樣
+    REFERENCE_FROM      //!< 參照, 執行完兩個Self會完全等效
+  };
 private:
   class Body {
   private:
     struct Kernel {
       Data    data_;
       int  counter_;
-      Kernel() {
-        counter_ = 1;
-      }
-      Kernel(Data const& data) {
-        counter_ = 1;
-        data_.copyFrom(data);
-      }
+
+      Kernel(                ): data_(    ), counter_(1) { }
+      Kernel(Data const& data): data_(data), counter_(1) { }
+      ~Kernel() { }
     };
-    Kernel *pointer_;
+
+    Kernel* pointer_;
     int     counter_;
   public:
-    Body() {
-      counter_ = 1;
-      pointer_ = new Kernel;
+    Body(             ): pointer_(new Kernel( )), counter_(1) { }
+    Body(Data const& d): pointer_(new Kernel(d)), counter_(1) { }
+    Body(Body const& b): pointer_(b.pointer_   ), counter_(1) {
+      ++pointer_->counter_;
+    }
+    ~Body() {
+      clear();
     }
-    Body(Body const& b) {
-      counter_ = 1;
+    Body& copyFrom(Body const& b) {
+      clear();
       pointer_ = b.pointer_;
-      pointer_->counter_++;
+      ++(pointer_->counter_);
+      return *this;
     }
-    ~Body() {
-      pointer_->counter_--;
+    void clear() {
+      --(pointer_->counter_);
       if (pointer_->counter_ <= 0) {
         delete pointer_;
       }
     }
-    int attatch() { return ++counter_; }
-    int detatch() { return --counter_; }
-    Data const* access() const { return &(pointer_->data_); }
-    Data      * modify() {
+    int attach() { return ++counter_; }
+    int detach() { return --counter_; }
+    Data const* access() const {
+      return &(pointer_->data_);
+    }
+    Data* modify() {
       if (pointer_->counter_ > 1) {
-        pointer_->counter_--;
+        --pointer_->counter_;
         pointer_ = new Kernel(pointer_->data_);
       }
       return &(pointer_->data_);
     }
   };
-  Body*  body_;
-  
-  void clear(Body* body) {
-    if (body != NULL) {
-      if (body->detatch() <= 0) {
-        delete body;
-      }
+
+  Body* body_;
+
+  void clear() {
+    if (body_->detach() <= 0) {
+      delete body_;
     }
   }
 public:
   /*!
-   *@brief constructor
+   *@brief constructor, 並宣告一個實體
+   */
+  Self(): body_(new Body()) {
+  }
+
+  /*!
+   * @brief connstructor, 宣告一個實體, 其中該實體的constructor用copy
+   *        constructor
+   *
+   * @param [in] d Inital data
    *
-   *@param [in] create_body 是否要new一個實體資料 (如果constructor完, 馬上就要
-   *                        \c copyFrom() , 或是 \c referenceFrom() 的話
-   *                        不太需要 new一個實體, 否則基本上都要
    */
-  Self(bool create_body) {
-    body_ = (create_body ? new Body() : NULL);
+  Self(Data const& d): body_(new Body(d)) {
+  }
+
+  /*!
+   * @brief constructor, 使用給定的Self當作init value, 並且可以指定要用
+   *         reference還是copy
+   *
+   * @param [in] b 給定的Self
+   * @param [in] d 指定要用copy還是reference
+   */
+  Self(Self const& b, DuplicateType d) {
+    switch(d) {
+    case COPY_FROM:
+      body_ = new Body(*b.body_);
+      break;
+    case REFERENCE_FROM:
+      body_ = b.body_;
+      body_->attach();
+      break;
+    }
   }
-  
-  //! @brief 不允許copy constructor
+
+  //! @brief Disallow copy constructor
   Self(Self const& b);
-  
+
   //! @brief 解構子
   ~Self() {
-    clear(body_);
+    clear();
   }
-  
+
   //! @brief 回傳指向 Data const 的指標
   Data const* operator->() const {
     return body_->access();
   }
-  
+
   //! @brief 回傳指向 Data 的指標, 如有需要, 這邊會做資料的duplicate
   Data* operator->() {
     return body_->modify();
   }
-  
+
   //! @brief 回傳非const型態的自己
   Self& operator()() const {
     return *((Self*)this);
   }
-  
+
   /*!
-   *@brief 將給定的 \c Self 的資料複製到自己這裡
+   * @brief 將給定的 \c Self 的資料複製到自己這裡
    *
-   *@param [in] s 給定的\c Self
-   *@return 無
+   * @param [in] s 給定的\c Self
+   * @return *this
    *
-   *@note 與reference的差別是, copy之後若該給定的 \c Self 有資料修改, 
-   *      this 這邊 \b 不會 被改到
+   * @note 與reference的差別是, copy之後若該給定的 \c Self 有資料修改,
+   *       this 這邊 \b 不會 被改到
    */
-  void copyFrom(Self const& s) {
-    Body* old = body_;
-    body_ = new Body(*(s.body_));
-    clear(old);
+  Self const& copyFrom(Self const& s) {
+    if (body_->access() != s.body_->access()) {
+      body_->copyFrom(*s.body_);
+    }
+    return *this;
   }
-  
+
   /*!
-   *@brief 將自己reference 到給定的 \c Self
+   * @brief 將自己reference 到給定的 \c Self
    *
-   *@param [in] s 給定的\c Self
-   *@return 無
+   * @param [in] s 給定的\c Self
+   * @return *this
    *
-   *@note 把reference想像成指標會比較容易思考, 譬如 \c a.referenceFrom(b)
-   *      \c b.referenceFrom(c) 相當於 \b a指向b原本指的記憶體位置,
-   *      \b b指向c原本指的記憶體位置 , 之後更動c時, 只有b會被牽連
-   *      
+   * @note 把reference想像成指標會比較容易思考, 譬如 \c a.referenceFrom(b)
+   *       \c b.referenceFrom(c) 相當於 \b a指向b原本指的記憶體位置,
+   *       \b b指向c原本指的記憶體位置 , 之後更動c時, 只有b會被牽連
    */
-  void referenceFrom(Self const& s) {
+  Self const& referenceFrom(Self const& s) {
     if (body_ != s.body_) {
-      clear(body_);
+      clear();
       body_ = s.body_;
-      body_->attatch();
+      body_->attach();
     }
+    return *this;
   }
-  
+
+  /*!
+   * @brief 依據給定DuplicateType選擇要呼叫copyFrom還是referenceFrom
+   *
+   * @param [in] s 給定的 \c Self
+   * @param [in] t 給定的 \c DuplicateType
+   * @return *this
+   */
+  Self const& duplicateFrom(Self const& s, DuplicateType t) {
+    switch(t) {
+    case COPY_FROM     : return      copyFrom(s);
+    case REFERENCE_FROM: return referenceFrom(s);
+    }
+    return *this;
+  }
+
   /*!
    * @brief 比對兩個 \c Self 是否指向同一個reference
    *
@@ -197,10 +239,10 @@ public:
   bool same(Self const& s) const {
     return (body_ == s.body_);
   }
-  
+
   /*!
    * @brief 比對兩個 \c Self 的內容是否一樣
-   * 
+   *
    * @param [in] s 另一個 \c Self
    * @return \c true/false 表示兩個內容是否一樣
    *
@@ -208,9 +250,9 @@ public:
    */
   bool equal(Self const& s) const {
     if (same(s) || body_->access() == s.body_->access()) return true;
-    return (body_->access()->equal(*(s.body_->access())));
+    return (*body_->access() == *s.body_->access());
   }
-  
+
   /*!
    * @brief 以reference作為判斷依據的小於判斷
    *
@@ -220,7 +262,7 @@ public:
   bool referenceLess(Self const& s) const {
     return (body_ < s.body_);
   }
-  
+
   //! @brief 將 \c operator= 給disable掉
   void operator=(Self const& a);
 };
diff --git a/meowpp/Usage.h b/meowpp/Usage.h
index 4537202..f0d2fc5 100644
--- a/meowpp/Usage.h
+++ b/meowpp/Usage.h
@@ -5,18 +5,14 @@
 
 #include <cstdlib>
 
+#include <algorithm>
 #include <string>
 #include <vector>
 #include <map>
-#include <algorithm>
-
-extern "C" {
-#include <unistd.h>
-}
 
 namespace meow {
 /*!
- *  @brief 管理參數設置, 自訂usage document, 分析argc, argv 
+ *  @brief 管理參數設置, 自訂usage document, 分析argc, argv
  *
  *  \b Usage 是用來分析argc, argv和輸出usage document的class. \n
  *  argc, argv的部份, 有以下規則
@@ -34,11 +30,10 @@ private:
   //! 存 (value, description)
   class Value {
   private:
-    String value_;
+    String       value_;
     String description_;
   public:
-    Value() {
-    }
+    Value() { }
     Value(String const& value, String const& description) {
       value_       = value;
       description_ = stringReplace(description, "<value>", value);
@@ -70,12 +65,12 @@ private:
     Option() {
     }
     Option(String const& description) {
-      has_setup_    = false;
-      has_value_    = false;
+      has_setup_   = false;
+      has_value_   = false;
       description_ = description;
-      must_setup_   = false;
+      must_setup_  = false;
     }
-    Option(String const& description, 
+    Option(String const& description,
            String const& type,
            String const& default_value,
            bool          must) {
@@ -110,31 +105,32 @@ private:
     bool valueAcceptAdd(String const& value, String const& description) {
       if (!has_value_) return false;
       if (std::find(values_accept_.begin(), values_accept_.end(),
-                    Value(value, "")) == values_accept_.end()){
+                    Value(value, "")) == values_accept_.end()) {
         values_accept_.push_back(Value(value, description));
       }
       return true;
     }
-    bool valueAcceptChk(String const& value){
+    bool valueAcceptChk(String const& value) {
       if (!has_value_) return false;
-      if (values_accept_.size() == 0) return true;
-      return (std::find(values_accept_.begin(), values_accept_.end(),
+      return (values_accept_.empty() ||
+              std::find(values_accept_.begin(), values_accept_.end(),
                         Value(value, "")) != values_accept_.end());
     }
     bool hasSetup() const{ return has_setup_; }
     bool hasValue() const{ return has_value_; }
     bool chkSetup() const{ return !(must_setup_ && !has_setup_); }
 
-    String usage(unsigned char opt, bool detail) const {
-      String ret(stringPrintf("-%c ", opt));
+    String usage(String opt, bool detail) const {
+      String ret(stringPrintf("%s%s ", (opt.size()>1 ? "--":"-"), opt.c_str()));
       if (!detail) {
         if (has_value_) ret += value_type_;
         if (!must_setup_) ret = "[" + ret + "]";
-      } else {
+      }
+      else {
         if (has_value_) {
           ret += value_type_ + " ";
           String default_string("");
-          if (value_default_ != "") 
+          if (value_default_ != "")
             default_string = "defalut='" + value_default_ + "'";
           String optional_string("");
           if (!must_setup_)
@@ -143,7 +139,8 @@ private:
           if (default_string.size() + optional_string.size() > 0) {
             if (default_string.size() > 0 && optional_string.size() > 0) {
               ret += "(" + optional_string + ", " + default_string + ")";
-            } else {
+            }
+            else {
               ret += "(" + optional_string + default_string + ")";
             }
           }
@@ -152,8 +149,7 @@ private:
         String accept_string;
         for (size_t i = 0; i < values_accept_.size(); i++) {
           if (i > 0)
-            accept_string += (i + 1 < values_accept_.size()
-                              ? ", " : " or ");
+            accept_string += (i + 1 < values_accept_.size() ? ", " : " or ");
           accept_string += "'" + values_accept_[i].value() + "'";
         }
         if (accept_string.size() == 0) accept_string = "... (anything)";
@@ -170,12 +166,14 @@ private:
       return ret;
     }
   };
-  typedef std::map<unsigned char, Option> Options;
-  typedef Options::const_iterator OptionsIterator;
-  String  name_;
-  Options options_;
-  Strings usage_begin_;
-  Strings usage_end_;
+  typedef std::map<String, Option> Options;
+  typedef Options::      iterator  OptionsIterator;
+  typedef Options::const_iterator  OptionsIteratorK;
+
+  String            name_;
+  Options        options_;
+  Strings    usage_begin_;
+  Strings      usage_end_;
   Strings proc_arguments_;
 public:
   /*!
@@ -216,14 +214,14 @@ public:
    *  @return \c true/false 表示 \b 是否成功
    */
   bool import(Usage const& usage) {
-    for (OptionsIterator
+    for (OptionsIteratorK
          it = usage.options_.begin(); it != usage.options_.end(); ++it) {
       if (options_.find(it->first) != options_.end())
         return false;
     }
-    for (OptionsIterator
+    for (OptionsIteratorK
          it = usage.options_.begin(); it != usage.options_.end(); ++it) {
-      options_[it->first] = it->second;
+      options_.insert(std::pair<String, Option>(it->first, it->second));
     }
     for (size_t i = 0; i < usage.usage_begin_.size(); ++i)
       usage_begin_.push_back(usage.usage_begin_[i]);
@@ -239,11 +237,12 @@ public:
    *  @return \c true/false 表 \b 是否成功
    */
   bool update(Usage const& usage) {
-    for (OptionsIterator
+    for (OptionsIteratorK
          it = usage.options_.begin(); it != usage.options_.end(); ++it) {
-      if (options_.find(it->first) == options_.end()) continue;
-      for(size_t i = 0, I = it->second.values().size(); i < I; i++){
-        options_[it->first].valueAdd(it->second.value(i));
+      OptionsIterator my = options_.find(it->first);
+      if (my == options_.end()) continue;
+      for (size_t i = 0, I = it->second.values().size(); i < I; ++i) {
+        my->second.valueAdd(it->second.value(i));
       }
     }
     return true;
@@ -256,9 +255,9 @@ public:
    *  @param [in] des 即description, 用來解釋這個選項的意義用的
    *  @return \c true/false 表 \b 是否成功
    */
-  bool optionAdd(unsigned char opt, String const& des) {
+  bool optionAdd(String opt, String const& des) {
     if (options_.find(opt) != options_.end()) return false;
-    options_[opt] = Option(des);
+    options_.insert(std::pair<String, Option>(opt, Option(des)));
     return true;
   }
 
@@ -272,12 +271,13 @@ public:
    *  @param [in] must 表示是否一定要設定
    *  @return \c true/false 表 \b 是否成功
    */
-  bool optionAdd(unsigned char opt, String const& des,
+  bool optionAdd(String opt, String const& des,
                  String const& val_type,
                  String const& val_default,
                  bool must) {
     if (options_.find(opt) != options_.end()) return false;
-    options_[opt] = Option(des, val_type, val_default, must);
+    options_.insert(std::pair<String, Option>(
+        opt, Option(des, val_type, val_default, must)));
     return true;
   }
 
@@ -289,11 +289,10 @@ public:
    *  @param [in] des 關於此額外選項的說明
    *  @return \c true/false 表 \b 是否成功
    */
-  bool optionValueAcceptAdd(unsigned char opt,
-                            String const& val,
-                            String const& des) {
-    if (options_.find(opt) == options_.end()) return false;
-    return options_[opt].valueAcceptAdd(val, des);
+  bool optionValueAcceptAdd(String opt, String const& val, String const& des) {
+    OptionsIterator o = options_.find(opt);
+    if (o == options_.end()) return false;
+    return o->second.valueAcceptAdd(val, des);
   }
 
   /*!
@@ -302,9 +301,9 @@ public:
    *  @param [in] opt 指定字符
    *  @return \c true/false 表 \b 是否有設定此選項
    */
-  bool hasOptionSetup(unsigned char opt) const {
-    return (options_.find(opt) != options_.end() &&
-            options_.find(opt)->second.hasSetup());
+  bool hasOptionSetup(String opt) const {
+    OptionsIteratorK o = options_.find(opt);
+    return (o != options_.end() && o->second.hasSetup());
   }
 
   /*!
@@ -313,9 +312,10 @@ public:
    *  @param [in] opt 指定字符
    *  @return 回傳次數
    */
-  size_t optionValuesSize(unsigned char opt) const {
-    if(options_.find(opt) == options_.end()) return 0;
-    return options_.find(opt)->second.values().size();
+  size_t optionValuesSize(String opt) const {
+    OptionsIteratorK o = options_.find(opt);
+    if (o == options_.end()) return 0;
+    return o->second.values().size();
   }
 
   /*!
@@ -325,11 +325,10 @@ public:
    *  @param [in] index 第幾個
    *  @return 回傳參數 \b -(opt) 的第 \b index 個額外選項
    */
-  String optionValue(unsigned char opt, size_t index) const {
-    if (options_.find(opt) == options_.end()) {
-      return String();
-    }
-    return options_.find(opt)->second.value(index);
+  String optionValue(String opt, size_t index) const {
+    OptionsIteratorK o = options_.find(opt);
+    if (o == options_.end()) return String();
+    return o->second.value(index);
   }
 
   /*!
@@ -386,16 +385,16 @@ public:
    *
    *  @return \b usage \b string
    */
-  String usage() const{
+  String usage() const {
     Usage::String out = stringPrintf("USAGE\n    %s", name_.c_str());
-    for (OptionsIterator
+    for (OptionsIteratorK
          it = options_.begin(); it != options_.end(); ++it)
       out += " " + it->second.usage(it->first, false);
     out += "\n\nDESCRIPTION\n";
     for (size_t i = 0; i < usage_begin_.size(); ++i) {
       out += "    " + usage_begin_[i] + "\n\n";
     }
-    for (OptionsIterator
+    for (OptionsIteratorK
          it = options_.begin(); it != options_.end(); ++it) {
       out += it->second.usage(it->first, true);
     }
@@ -407,52 +406,65 @@ public:
 
   /*!
    *  @brief 給定argc, argv, 將各參數設置
-   *  @param [in] argc,argv 
+   *  @param [in]  argc,argv
    *  @param [out] errmsg 將錯誤訊息寫到這裡
    *                     (若給定NULL pointer, 則會把錯誤訊息忽略)
    *  @return \c true/false \b 成功與否 (否的話代表有錯誤的設定值在其中)
    */
-  bool arguments(int argc, char** argv, String* errmsg){
-    opterr = 0;
-    String s;
-    OptionsIterator it;
-    String zzz;
+  bool arguments(int argc, char** argv, String* errmsg) {
+    String  zzz;
     String& err = (errmsg == NULL ? zzz : *errmsg);
-    for (it = options_.begin(); it != options_.end(); ++it) {
-      s += (char)(it->first);
-      if (it->second.hasValue()) s += ":";
-    }
-    bool succ = true;
-    for (int opt; (opt = getopt(argc, argv, s.c_str())) != -1; ) {
-      if (options_.find(opt) == options_.end()) {
-        if(options_.find(optopt) == options_.end()){
-          err += stringPrintf("Unknown option '-%c'\n", optopt);
-        }else{
-          err += stringPrintf("No specify argument to '-%c'\n",
-                              optopt);
+    for (int i = 0; i < argc; ++i) {
+      if (argv[i][0] == '-') {
+        String opt;
+        if (argv[i][1] == '-') { // long option
+          if (argv[i][2] == '\0' || argv[i][3] == '\0') {
+            err = stringPrintf("Invalid option '%s'", argv[i]);
+            return false;
+          }
+          opt = String(argv[i] + 2);
+        }
+        else { // short option
+          if (argv[i][1] == '\0' || argv[i][2] != '\0') {
+            err = stringPrintf("Invalid option '%s'", argv[i]);
+            return false;
+          }
+          opt = String(argv[i] + 1);
+        }
+        OptionsIterator o = options_.find(opt);
+        if (o == options_.end()) {
+          err = stringPrintf("Unknown option '%s'", argv[i]);
+          return false;
+        }
+        if (o->second.hasValue()) {
+          if (i + 1 >= argc) {
+            err = stringPrintf("Option '%s' need a value", argv[i]);
+            return false;
+          }
+          if (o->second.valueAdd(argv[i + 1]) < 0) {
+            err = stringPrintf("Invalid value '%s' of options '%s'",
+                               argv[i + 1], argv[i]);
+            return false;
+          }
+          ++i;
+        }
+        else {
+          o->second.valueAdd("");
         }
-        succ = false;
-        continue;
       }
-      if (options_[opt].valueAdd(optarg == NULL ? "" : optarg) < 0) {
-        err += stringPrintf("Option argument '%s' to '-%c' is not allowed\n"
-                            , optarg, opt);
-        succ = false;
-        continue;
+      else {
+        proc_arguments_.push_back(String(argv[i] + (argv[i][0]=='\\' ? 1:0)));
       }
     }
-    for (it = options_.begin(); it != options_.end(); it++) {
+    for (OptionsIteratorK it = options_.begin(); it != options_.end(); ++it) {
       if (it->second.chkSetup() == false) {
-        err += stringPrintf("No specify argument to '-%c'\n",
-                            it->first);
-        succ = false;
-        continue;
+        err += stringPrintf("No specify argument to '%s%s'\n",
+                            (it->first.size() > 1 ? "--" : "-"),
+                            it->first.c_str());
+        return false;
       }
     }
-    for (int i = optind; i < argc; i++) {
-      proc_arguments_.push_back(String(argv[i]));
-    }
-    return succ;
+    return true;
   }
 };
 
diff --git a/meowpp/colors/Color3_Space.h b/meowpp/colors/Color3_Space.h
index 87e215b..94cae39 100644
--- a/meowpp/colors/Color3_Space.h
+++ b/meowpp/colors/Color3_Space.h
@@ -1,18 +1,17 @@
 #ifndef   colors_Color3_Space_H__
 #define   colors_Color3_Space_H__
 
-#include <cstdlib>
-
-#include "../Self.h"
 #include "../geo/Vectors.h"
 #include "../math/Matrix.h"
 #include "../math/utility.h"
 
+#include <cstdlib>
+
 namespace meow {
 
 /*!
  * @brief 以三個channel所組成的色彩空間
- * 
+ *
  * @author cat_leopard
  */
 template<class T>
@@ -31,21 +30,14 @@ protected:
   Color3_Space(Color3_Space const& b) {
     min_ = b.min_;
     max_ = b.max_;
-    copyFrom(b);
+    val_ = b.val_;
   }
-public:
-  virtual ~Color3_Space() { }
   Color3_Space<T>& copyFrom(Color3_Space<T> const& b) {
     val_ = b.val_;
     return *this;
   }
-  template<class T2>
-  Color3_Space<T>& copyFrom(Color3_Space<T2> const& b) {
-    for (size_t i = 0; i < 3; i++) {
-      val(i, ratioMapping((double)b.min(i), (double)b.max(i), (double)b.rgb(i),
-                          (double)min(i), (double)max(i)));
-    }
-  }
+public:
+  virtual ~Color3_Space() { }
   Vector3D<T> const& minV() const { return min_; }
   Vector3D<T> const& maxV() const { return max_; }
   Vector3D<T> const& valV() const { return val_; }
diff --git a/meowpp/colors/HSL_Space.h b/meowpp/colors/HSL_Space.h
index f70eab1..52c218f 100644
--- a/meowpp/colors/HSL_Space.h
+++ b/meowpp/colors/HSL_Space.h
@@ -2,10 +2,10 @@
 #define   colors_HSL_Space_H__
 
 #include "Color3_Space.h"
-#include "../geo/Vectors.h"
-
 #include "RGB_Space.h"
 #include "YUV_Space.h"
+
+#include "../geo/Vectors.h"
 #include "../math/utility.h"
 
 #include <cstdlib>
@@ -16,7 +16,7 @@ namespace meow {
  * @brief 以浮點數\b Y(亮度), \b U(色度), \b V(濃度) 三個值所組成的色彩空間
  *
  * 其中範圍都介於0.0~1.0之間
- * 
+ *
  * @author cat_leopard
  */
 class HSLf_Space: public Color3_Space<double> {
diff --git a/meowpp/colors/HSV_Space.h b/meowpp/colors/HSV_Space.h
index 96ba33c..eb8971b 100644
--- a/meowpp/colors/HSV_Space.h
+++ b/meowpp/colors/HSV_Space.h
@@ -17,7 +17,7 @@ namespace meow {
  * @brief 以浮點數\b Y(亮度), \b U(色度), \b V(濃度) 三個值所組成的色彩空間
  *
  * 其中範圍都介於0.0~1.0之間
- * 
+ *
  * @author cat_leopard
  */
 class HSVf_Space: public Color3_Space<double> {
diff --git a/meowpp/colors/RGB_Space.h b/meowpp/colors/RGB_Space.h
index ea04598..2523d8b 100644
--- a/meowpp/colors/RGB_Space.h
+++ b/meowpp/colors/RGB_Space.h
@@ -13,7 +13,7 @@ namespace meow {
  * @brief 以整數 \b Red, \b Green, \b Blue 三個值所組成的色彩空間
  *
  * 其中範圍都介於0~255之間
- * 
+ *
  * @author cat_leopard
  */
 class RGBi_Space: public Color3_Space<int> {
@@ -80,7 +80,7 @@ public:
  * @brief 以浮點數\b Red, \b Green, \b Blue 三個值所組成的色彩空間
  *
  * 其中範圍都介於0.0~1.0之間
- * 
+ *
  * @author cat_leopard
  */
 class RGBf_Space: public Color3_Space<double> {
diff --git a/meowpp/colors/YUV_Space.h b/meowpp/colors/YUV_Space.h
index 651caaf..8390034 100644
--- a/meowpp/colors/YUV_Space.h
+++ b/meowpp/colors/YUV_Space.h
@@ -15,7 +15,7 @@ namespace meow {
  * @brief 以浮點數\b Y(亮度), \b U(色度), \b V(濃度) 三個值所組成的色彩空間
  *
  * 其中範圍都介於0.0~1.0之間
- * 
+ *
  * @author cat_leopard
  */
 class YUVf_Space: public Color3_Space<double> {
diff --git a/meowpp/dsa/DisjointSet.h b/meowpp/dsa/DisjointSet.h
index 4575835..9be9c35 100644
--- a/meowpp/dsa/DisjointSet.h
+++ b/meowpp/dsa/DisjointSet.h
@@ -13,11 +13,11 @@ namespace meow {
  * 相關資料可參考
  * <a href="http://www.csie.ntnu.edu.tw/~u91029/DisjointSets.html">
  *   演算法筆記
- * </a> 
+ * </a>
  *
  * @note
  *    - 時間複雜度 \b 非常快 表示它真的算的超級快, 可以視為常數時間
- *    - 預設值所有 \a number 所在的集合的編號就是 \a number 本身, 
+ *    - 預設值所有 \a number 所在的集合的編號就是 \a number 本身,
  *      即沒有任兩個數在同一個set裡面
  *
  * @author cat_leopard
@@ -119,7 +119,7 @@ public:
    *
    * 將 \a number1 所在的集合 跟 \b number2 所在的集合 \b 合併,
    * 並回傳合併後新的集合的編號. \n
-   * 時間複雜度\b 非常快 
+   * 時間複雜度\b 非常快
    *
    * @param [in] a 即上述\a number1
    * @param [in] b 即上述\a number2
diff --git a/meowpp/dsa/HashTable.h b/meowpp/dsa/HashTable.h
index 9171c72..5f343f5 100644
--- a/meowpp/dsa/HashTable.h
+++ b/meowpp/dsa/HashTable.h
@@ -25,7 +25,7 @@ public:
 
   /*!
    * @brief constructor
-   * 
+   *
    * 設定table size, hash function
    */
   HashTableList(size_t size, HashFunc const& func): table_(size), func_(func) {
@@ -193,18 +193,18 @@ public:
     }
     return ret;
   }
-  
+
   //! @brief same as \c copyFrom(h)
   HashTableList& operator=(HashTableList const& h) {
     return copyFrom(h);
   }
-  
+
   //! @brief same as \c add(h)
   HashTableList& operator+=(HashTableList const& h) {
     add(h);
     return *this;
   }
-  
+
   //! @brief same as \c del(h)
   HashTableList& operator-=(HashTableList const& h) {
     del(h);
diff --git a/meowpp/dsa/KD_Tree.h b/meowpp/dsa/KD_Tree.h
index 05f9b1b..e3bd73b 100644
--- a/meowpp/dsa/KD_Tree.h
+++ b/meowpp/dsa/KD_Tree.h
@@ -43,12 +43,12 @@ private:
     Vector  vector_;
     ssize_t lChild_;
     ssize_t rChild_;
-    
+
     Node(Vector v, ssize_t l, ssize_t r): vector_(v), lChild_(l), rChild_(r){
     }
   };
   typedef std::vector<Node> Nodes;
-  
+
   class Sorter {
   private:
     Nodes const* nodes_;
@@ -187,16 +187,16 @@ private:
 public:
   //! Custom Type: Vectors is \c std::vector<Vector>
   typedef typename std::vector<Vector> Vectors;
-  
+
   //! @brief constructor, with dimension = 1
   KD_Tree(): kNIL_(-1), root_(kNIL_), needRebuild_(false), dimension_(1) {
   }
-  
+
   //! @brief constructor, given dimension
   KD_Tree(size_t dimension):
   kNIL_(-1), root_(kNIL_), needRebuild_(false), dimension_(dimension) {
   }
-  
+
   //! @brief destructor
   ~KD_Tree() {
   }
@@ -254,7 +254,7 @@ public:
   }
 
   /*!
-   * @brief 查找 
+   * @brief 查找
    *
    * 於set中找尋距離指定向量前 \c i 近的向量, 並依照由近而遠的順序排序.
    * 如果有兩個向量\c v1,v2 距離一樣, 且 \c cmp 為\c true , 則直接依照
diff --git a/meowpp/dsa/MergeableHeap.h b/meowpp/dsa/MergeableHeap.h
index af7ad75..0967edd 100644
--- a/meowpp/dsa/MergeableHeap.h
+++ b/meowpp/dsa/MergeableHeap.h
@@ -7,7 +7,7 @@
 namespace meow {
 
 /*!
- * @brief 
+ * @brief
  *
  * 一個用 \b 左偏樹 實作的 \c Maximum-Heap , 除了原本heap有的功能外,
  * 還支援 \c merge 功能
diff --git a/meowpp/dsa/SegmentTree.h b/meowpp/dsa/SegmentTree.h
index b2fa749..64eab4c 100644
--- a/meowpp/dsa/SegmentTree.h
+++ b/meowpp/dsa/SegmentTree.h
@@ -37,7 +37,7 @@ namespace meow {
  * - 若要維護區間最總和, 即每次都是詢問範圍 `[a, b]` 的總和, 則可以定義
  *      - \c operator+ 為 '回傳相加值'
  *      - \c operator* 為 '回傳(*this) * n'
- *      - \c operator| 為 '回傳相加值' 
+ *      - \c operator| 為 '回傳相加值'
  *
  * @author cat_leopard
  */
@@ -140,7 +140,7 @@ public:
     nodes_ = b.nodes_;
     return *this;
   }
-  
+
   /*!
    * @brief 回傳size
    */
@@ -182,7 +182,7 @@ public:
     if (rangeCorrect(&first, &last) == false) return ;
     update(first, last, 0, size_ - 1, 0, delta, false);
   }
-  
+
   //! @brief same as copyFrom(b)
   SegmentTree& operator=(SegmentTree const& b) {
     return copyFrom(b);
diff --git a/meowpp/dsa/SplayTree.h b/meowpp/dsa/SplayTree.h
index 5e9fb3c..40a2a0b 100644
--- a/meowpp/dsa/SplayTree.h
+++ b/meowpp/dsa/SplayTree.h
@@ -43,7 +43,7 @@ private:
     size_t        size_;
     Node*       parent_;
     Node*        child_[2];
-    
+
     Node(Key const& key, Value const& value):
     key_(key), keyOffset_(0), value_(value) {
       size_     = 1;
@@ -373,7 +373,7 @@ public:
 
   /*!
    * @brief 回傳一個指向NULL的Element,
-   * 
+   *
    * 以供 \c find ,\c order ,\c first ,\c last 等判斷是否有找到相對應的Element
    */
   Element end() const {
@@ -404,7 +404,7 @@ public:
 
   /*!
    * @brief 插入一組\c (Key ---> \c Value)
-   * 
+   *
    * 檢查是否已有Element的Key 為 \c key, 若有則回傳 \c false , 否則將
    * 一個 (Key -> Value) = (\c key -> \c value)的Element加入, 並回傳 \c true
    */
@@ -504,7 +504,7 @@ public:
   /*!
    * @brief 合併
    *
-   * 檢查是否自己中的 Key 都小於 \c tree2 中的Key, 或是完全相反, 
+   * 檢查是否自己中的 Key 都小於 \c tree2 中的Key, 或是完全相反,
    * 是的話把 \c tree2`中的 Element 都搬到自己這,
    * 同時清空 \c tree2 , 否則回傳 \c false
    */
@@ -522,7 +522,7 @@ public:
     tree2->root_ = NULL;
     return true;
   }
-  
+
   /*!
    * @brief 就像\c stl::map::operator[]
    *
@@ -578,7 +578,7 @@ private:
     size_t        size_;
     Node*       parent_;
     Node*        child_[2];
-    
+
     Node(Key const& key, Value const& value):
     valueOffset_(0), range_(value),
     key_(key), keyOffset_(0), value_(value) {
@@ -932,7 +932,7 @@ public:
 
   /*!
    * @brief 回傳一個指向NULL的Element,
-   * 
+   *
    * 以供 \c find ,\c order ,\c first ,\c last 等判斷是否有找到相對應的Element
    */
   Element end() const {
@@ -952,7 +952,7 @@ public:
   bool empty() const{
     return (size() == 0);
   }
-  
+
   /*!
    * @brief 查找
    *
@@ -992,7 +992,7 @@ public:
 
   /*!
    * @brief 插入一組\c (Key ---> \c Value)
-   * 
+   *
    * 檢查是否已有Element的Key 為 \c key, 若有則回傳 \c false , 否則將
    * 一個 (Key -> Value) = (\c key -> \c value)的Element加入, 並回傳 \c true
    */
@@ -1110,7 +1110,7 @@ public:
   /*!
    * @brief 合併
    *
-   * 檢查是否自己中的 Key 都小於 \c tree2 中的Key, 或是完全相反, 
+   * 檢查是否自己中的 Key 都小於 \c tree2 中的Key, 或是完全相反,
    * 是的話把 \c tree2`中的 Element 都搬到自己這,
    * 同時清空 \c tree2 , 否則回傳 \c false
    */
@@ -1128,7 +1128,7 @@ public:
     tree2->root_ = NULL;
     return true;
   }
-  
+
   /*!
    * @brief 就像\c stl::map::operator[]
    *
diff --git a/meowpp/dsa/VP_Tree.h b/meowpp/dsa/VP_Tree.h
index 75186e6..9c85930 100644
--- a/meowpp/dsa/VP_Tree.h
+++ b/meowpp/dsa/VP_Tree.h
@@ -15,13 +15,13 @@ namespace meow {
 /*!
  * @brief 跟KD_Tree很像歐
  *
- * \c VP_Tree 用來維護由 \b N個K維度向量所成的集合 , 
+ * \c VP_Tree 用來維護由 \b N個K維度向量所成的集合 ,
  * 並可於該set中查找 \b 前i個離給定向量最接近的向量* .
- * 不像 \c KD_Tree 二分樹每次都選擇一個維度去分, 分成小的跟大的, 
+ * 不像 \c KD_Tree 二分樹每次都選擇一個維度去分, 分成小的跟大的,
  * \c VP_Tree 每次選一個點, 將資料分成 離這個點近的, 跟離這個點遠的.
  * 至於怎麼選呢...., 嘛還沒研究, 先random
  *
- * 參考資料連結: 
+ * 參考資料連結:
  *   - http://stevehanov.ca/blog/index.php?id=130
  *   - http://pnylab.com/pny/papers/vptree/vptree
  *
@@ -33,7 +33,7 @@ namespace meow {
  * |const |  Vector|operator[] |(size_t \c n) | Scalar     | 取得第\c n 維度量 |
  * |const |  Vector|operator=  |(Vector \c v) | Vector&    | copy operator     |
  * |const |  Vector|operator<  |(Vector \c v) | bool       | 權重比較          |
- * |const |  Scalar| 'Scalar'  |(int    \c n) | Scalar     | 建構子,  
+ * |const |  Scalar| 'Scalar'  |(int    \c n) | Scalar     | 建構子,
  *                                                           其中一定\c n=0or4 |
  * |const |  Scalar|operator*  |(Scalar \c s) | Scalar     | 相乘              |
  * |const |  Scalar|operator+  |(Scalar \c s) | Scalar     | 相加              |
@@ -87,12 +87,12 @@ private:
   };
   typedef std::vector<Answer> AnswerV;
   typedef std::priority_queue<Answer, AnswerV, AnswerCompare> Answers;
-  
+
   Vectors     vectors_;
   Node*          root_;
   size_t    dimension_;
   bool    needRebuild_;
-  
+
   Scalar distance2(Vector const& v1, Vector const& v2) const {
     Scalar ret(0);
     for (size_t i = 0; i < dimension_; i++) ret += squ(v1[i] - v2[i]);
@@ -212,7 +212,7 @@ public:
   VP_Tree(): root_(NULL), vectors_(0), dimension_(1), needRebuild_(false){
     reset(0);
   }
-  
+
   //! @brief constructor, 複製資料
   VP_Tree(VP_Tree const& tree2):
   vectors_(tree2.vectors_),
@@ -220,7 +220,7 @@ public:
   dimension_(tree2.dimension_),
   needRebuild_(tree2.needRebuild_) {
   }
-  
+
   //! @brief constructor, 給定dimension
   VP_Tree(size_t dimension):
   vectors_(0),
@@ -229,12 +229,12 @@ public:
   needRebuild_(false) {
     reset(dimension);
   }
-  
+
   //! @brief destructor
   ~VP_Tree() {
     clear(root_);
   }
-  
+
   /*!
    * @brief 複製資料
    */
@@ -287,7 +287,7 @@ public:
   }
 
   /*!
-   * @brief 查找 
+   * @brief 查找
    *
    * 於set中找尋距離指定向量前 \c i 近的向量, 並依照由近而遠的順序排序.
    * 如果有兩個向量\c v1,v2 距離一樣, 且 \c cmp 為\c true , 則直接依照
@@ -325,7 +325,7 @@ public:
     dimension_ = std::max((size_t)1, dimension);
     return dimension_;
   }
-  
+
   //! @brief same as \c copyFrom(tree2)
   VP_Tree& operator=(VP_Tree const& tree2) {
     return copyFrom(tree2);
diff --git a/meowpp/geo/Vectors.h b/meowpp/geo/Vectors.h
index 8cdafc9..0f09ee6 100644
--- a/meowpp/geo/Vectors.h
+++ b/meowpp/geo/Vectors.h
@@ -22,165 +22,165 @@ public:
   //! @brief consturctor (0, 0)
   Vector2D(): x_(0), y_(0) {
   }
-  
+
   //! @brief consturctor (from another Vector2D)
   Vector2D(Vector2D const& v): x_(v.x_), y_(v.y_) {
   }
-  
+
   //! @brief constructor (s, s)
   Vector2D(Scalar const& s): x_(s), y_(s) {
   }
-  
+
   //! @brief constructor (sx, sy)
   Vector2D(Scalar const& sx, Scalar const& sy): x_(sx), y_(sy) {
   }
-  
+
   //! @brief constructor (from another Vector)
   Vector2D(Vector<Scalar> const& v): x_(v(0)), y_(v(1)) {
   }
-  
+
   //! @brief constructor (from another Vector, i-th)
   Vector2D(Vector<Scalar> const& v, size_t i): x_(v(i)), y_(v(i + 1)) {
   }
-  
+
   //! @brief destructor
   ~Vector2D() {
   }
-  
+
   //! @brief copy
   Vector2D& copyFrom(Vector2D const& v) {
     return xy(v.x(), v.y());
   }
-  
+
   //! @brief access x
   Scalar const& x() const {
     return x_;
   }
-  
+
   //! @brief access x with non constant reference
   Scalar& xGet() {
     return x_;
   }
-  
+
   //! @brief access y with non constant reference
   Scalar& yGet() {
     return y_;
   }
-  
+
   //! @brief access y
   Scalar const& y() const {
     return y_;
   }
-  
+
   //! @brief modify x
   Scalar const& x(Scalar const& s) {
     x_ = s;
     return x();
   }
-  
+
   //! @brief modify y
   Scalar const& y(Scalar const& s) {
     y_ = s;
     return y();
   }
-  
+
   //! @brief modify x and y
   Vector2D& xy(Scalar const& sx, Scalar const& sy){
     x(sx);
     y(sy);
     return *this;
   }
-  
+
   //! @brief return \a +(*this)
   Vector2D positive() const {
     return *this;
   }
-  
+
   //! @brief return \a -(*this)
   Vector2D negative() const {
     return Vector2D(-x(), -y());
   }
-  
+
   //! @brief return \a count-clockwise \a rotate \a 90 \a degree of itself
   Vector2D right()const{
     return Vector2D(-y(), x());
   }
-  
+
   //! @brief return \a (*this)+v
   Vector2D add(Vector2D const& v) const {
     return Vector2D(x() + v.x(), y() + v.y());
   }
-  
+
   //! @brief Let itself add v
   Vector2D& added(Vector2D const& v) {
     return xy(x() + v.x(), y() + v.y());
   }
-  
+
   //! @brief return \a (*this)-v
   Vector2D sub(Vector2D const& v) const {
     return Vector2D(x() - v.x(), y() - v.y());
   }
-  
+
   //! @brief Let itself substract v
   Vector2D& subed(Vector2D const& v) {
     return xy(x() - v.x(), y() - v.y());
   }
-  
+
   //! @brief return \a (*this)*s , where s is a scalar
   Vector2D mul(Scalar const& s) const {
     return Vector2D(x() * s, y() * s);
   }
-  
+
   //! @brief Let itself mulitple s
   Vector2D& muled(Scalar const& s) {
     return xy(x() * s, y() * s);
   }
-  
+
   //! @brief return \a (*this)/s , where s is a scalar
   Vector2D div(Scalar const& s) const {
     return Vector2D(x() / s, y() / s);
   }
-  
+
   //! @brief Let itself divide s
   Vector2D& dived(Scalar const& s) {
     return xy(x() / s, y() / s);
   }
-  
+
   //! @brief same as dot(v)
   Scalar mul(Vector2D const& v) const {
     return dot(v);
   }
-  
+
   //! @brief dot
   Scalar dot(Vector2D const& v) const {
     return x() * v.x() + y() * v.y();
   }
-  
+
   //! @brief cross
   Scalar cross(Vector2D const& v) const {
     return x() * v.y() - y() * v.x();
   }
-  
+
   //! @brief sqrt of length2
   Scalar length() const {
     return Scalar(sqrt(double(length2())));
   }
-  
+
   //! @brief same as \a dot(*this)
   Scalar length2() const {
     return dot(*this);
   }
-  
+
   //! @brief return normalize form of itself
   Vector2D normalize() const {
     return div(length());
   }
-  
+
   //! @brief normalize itself
   Vector2D& normalized() {
     return dived(length());
   }
-  
+
   //! @brief return rotate \a theta degree of itself
   Vector2D rotate(Scalar const& theta) const {
     Scalar cs(cos(-double(theta)));
@@ -188,22 +188,22 @@ public:
     Vector2D<Scalar> new_x(cs, sn);
     return Vector2D(new_x.dot(*this), new_x.cross(*this));
   }
-  
+
   //! @brief Let itself rotate \a theta degree
   Vector2D& rotated(Scalar const& theta) {
     return copyFrom(rotate(theta));
   }
-  
+
   //! @brief return reflect from given vector \a v
   Vector2D reflect(Vector2D const& v) const {
     return v.mul(v.dot(*this) * 2 / v.length2()).sub(*this);
   }
-  
+
   //! @brief reflect itself given vector \a v
   Vector2D& reflected(Vector2D const& v) {
     return copyFrom(reflecte(v));
   }
-  
+
   //! @brief return a 2x1 matrix form of itself
   Matrix<Scalar> matrix() const {
     static Matrix<Scalar> ret(2, 1, Scalar(0));
@@ -211,7 +211,7 @@ public:
     ret(1, 0, y());
     return ret;
   }
-  
+
   //! @brief return a 3x1 matrix form of itself
   Matrix<Scalar> matrix(Scalar const& homo) const {
     static Matrix<Scalar> ret(3, 1, Scalar(0));
@@ -220,11 +220,11 @@ public:
     ret(2, 0, homo);
     return ret;
   }
-  
+
   Scalar const& operator()(size_t n) const {
     return (n == 0 ? x() : y());
   }
-  
+
   Vector2D& operator()(Scalar const& sx, Scalar const& sy) {
     return xy(sx, sy);
   }
@@ -259,86 +259,86 @@ public:
   //! @brief consturctor (0, 0)
   Vector3D(): x_(0), y_(0), z_(0) {
   }
-  
+
   //! @brief consturctor (from another Vector3D)
   Vector3D(Vector3D const& v): x_(v.x_), y_(v.y_), z_(v.z_) {
   }
-  
+
   //! @brief constructor (s, s)
   Vector3D(Scalar const& s): x_(s), y_(s), z_(s) {
   }
-  
+
   //! @brief constructor (sx, sy)
   Vector3D(Scalar const& sx,
            Scalar const& sy,
            Scalar const& sz): x_(sx), y_(sy), z_(sz) {
   }
-  
+
   //! @brief constructor (from another Vector)
   Vector3D(Vector<Scalar> const& v): x_(v(0)), y_(v(1)), z_(v(2)) {
   }
-  
+
   //! @brief constructor (from another Vector, i-th)
   Vector3D(Vector<Scalar> const& v, size_t i): x_(v(i)), y_(v(i+1)), z_(v(i+2)){
   }
-  
+
   //! @brief destructor
   ~Vector3D(){
   }
-  
+
   //! @brief copy
   Vector3D& copyFrom(Vector3D const& v){
     return xyz(v.x(), v.y(), v.z());
   }
-  
+
   //! @brief access x
   Scalar const& x() const{
     return x_;
   }
-  
+
   //! @brief access y
   Scalar const& y() const{
     return y_;
   }
-  
+
   //! @brief access z
   Scalar const& z() const{
     return z_;
   }
-  
+
   //! @brief access x with non constant reference
   Scalar& xGet() {
     return x_;
   }
-  
+
   //! @brief access y with non constant reference
   Scalar& yGet() {
     return y_;
   }
-  
+
   //! @brief access z with non constant reference
   Scalar& zGet() {
     return z_;
   }
-  
+
   //! @brief modify x
   Scalar const& x(Scalar const& s) {
     x_ = s;
     return x();
   }
-  
+
   //! @brief modify y
   Scalar const& y(Scalar const& s) {
     y_ = s;
     return y();
   }
-  
+
   //! @brief modify z
   Scalar const& z(Scalar const& s) {
     z_ = s;
     return z();
   }
-  
+
   //! @brief modify x and y
   Vector3D& xyz(Scalar const& sx, Scalar const& sy, Scalar const& sz) {
     x(sx);
@@ -346,99 +346,99 @@ public:
     z(sz);
     return *this;
   }
-  
+
   //! @brief return \a +(*this)
   Vector3D positive() const {
     return *this;
   }
-  
+
   //! @brief return \a -(*this)
   Vector3D negative() const {
     return Vector3D(-x(), -y(), -z());
   }
-  
+
   //! @brief return \a (*this)+v
   Vector3D add(Vector3D const& v) const {
     return Vector3D(x() + v.x(), y() + v.y(), z() + v.z());
   }
-  
+
   //! @brief Let itself add v
   Vector3D& added(Vector3D const& v) {
     return xyz(x() + v.x(), y() + v.y(), z() + v.z());
   }
-  
+
   //! @brief return \a (*this)-v
   Vector3D sub(Vector3D const& v) const {
     return Vector3D(x() - v.x(), y() - v.y(), z() - v.z());
   }
-  
+
   //! @brief Let itself substract v
   Vector3D& subed(Vector3D const& v) {
     return xyz(x() - v.x(), y() - v.y(), z() - v.z());
   }
-  
+
   //! @brief return \a (*this)*s , where s is a scalar
   Vector3D mul(Scalar const& s) const {
     return Vector3D(x() * s, y() * s, z() * s);
   }
-  
+
   //! @brief Let itself mulitple s
   Vector3D& muled(Scalar const& s) {
     return xyz(x() * s, y() * s, z() * s);
   }
-  
+
   //! @brief return \a (*this)/s , where s is a scalar
   Vector3D div(Scalar const& s) const {
     return Vector3D(x() / s, y() / s, z() / s);
   }
-  
+
   //! @brief Let itself divide s
   Vector3D& dived(Scalar const& s) {
     return xyz(x() / s, y() / s, z() / s);
   }
-  
+
   //! @brief same as dot(v)
   Scalar mul(Vector3D const& v) const {
     return dot(v);
   }
-  
+
   //! @brief dot
   Scalar dot(Vector3D const& v) const {
     return x() * v.x() + y() * v.y() + z() * v.z();
   }
-  
+
   //! @brief cross
   Vector3D cross(Vector3D const& v) const {
     return Vector3D(y() * v.z() - z() * v.y(),
                     z() * v.x() - x() * v.z(),
                     x() * v.y() - y() * v.x());
   }
-  
+
   //! @brief crossed
   Vector3D& crossed(Vector3D const& v) {
     return copyFrom(cross(v));
   }
-  
+
   //! @brief sqrt of length2
   Scalar length() const {
     return Scalar(sqrt(double(length2())));
   }
-  
+
   //! @brief same as \a dot(*this)
   Scalar length2() const {
     return dot(*this);
   }
-  
+
   //! @brief return normalize form of itself
   Vector3D normalize() const {
     return div(length());
   }
-  
+
   //! @brief normalize itself
   Vector3D& normalized() {
     return dived(length());
   }
-  
+
   //! @brief return rotate \a theta degree by \a axis of itself
   Vector3D rotate(Vector3D const& axis, double theta) const {
     Vector3D a(axis.normalize());
@@ -446,22 +446,22 @@ public:
     Vector3D yy(a.cross(*this).mul(sin(theta)));
     return a.mul(a.dot(*this)).add(xx).add(yy);
   }
-  
+
   //! @brief Let itself rotate \a theta degree
   Vector3D& rotated(Vector3D const& axis, double theta) {
     return copyFrom(rotate(axis, theta));
   }
-  
+
   //! @brief return reflect from given vector \a v
   Vector3D reflect(Vector3D const& v) const {
     return v.mul(v.dot(*this) * 2 / v.length2()).sub(*this);
   }
-  
+
   //! @brief reflect itself given vector \a v
   Vector3D& reflected(Vector3D const& v) {
     return copyFrom(reflecte(v));
   }
-  
+
   //! @brief return a 3x1 matrix form of itself
   Matrix<Scalar> matrix() const {
     static Matrix<Scalar> ret(3, 1, Scalar(0));
@@ -470,7 +470,7 @@ public:
     ret(2, 0, z());
     return ret;
   }
-  
+
   //! @brief return a 3x1 matrix form of itself
   Matrix<Scalar> matrix(Scalar const& homo) const {
     static Matrix<Scalar> ret(4, 1, Scalar(0));
@@ -480,11 +480,11 @@ public:
     ret(3, 0, homo);
     return ret;
   }
-  
+
   Scalar const& operator()(size_t n) const {
     return (n == 0 ? x() : (n == 1 ? y() : z()));
   }
-  
+
   Vector3D& operator()(Scalar const& sx, Scalar const& sy, Scalar const& sz) {
     return xyz(sx, sy, sz);
   }
@@ -497,7 +497,7 @@ public:
   Vector3D operator*(Scalar   const& s) const { return mul(s); }
   Vector3D operator/(Scalar   const& s) const { return div(s); }
   Scalar   operator*(Vector3D const& v) const { return mul(v); }
-  
+
   Vector3D& operator=(Vector3D const& v) { return copyFrom(v); }
   Vector3D& operator+=(Vector3D const& v) { return added(v); }
   Vector3D& operator-=(Vector3D const& v) { return subed(v); }
diff --git a/meowpp/gra/Bitmap.h b/meowpp/gra/Bitmap.h
index c4ed4e0..59136c2 100644
--- a/meowpp/gra/Bitmap.h
+++ b/meowpp/gra/Bitmap.h
@@ -1,9 +1,6 @@
 #ifndef   gra_Bitmap_H__
 #define   gra_Bitmap_H__
 
-
-#include "../Self.h"
-
 #include "../math/utility.h"
 #include "../math/Matrix.h"
 
@@ -15,7 +12,7 @@
 #include <typeinfo>
 #include <cstdlib>
 
-namespace meow{
+namespace meow {
 
 /*!
  * @brief 二維點陣資料
@@ -26,7 +23,8 @@ template<class Pixel>
 class Bitmap: public ObjBase {
 private:
   Matrix<Pixel> matrix_;
-  //
+
+  //! 回傳高斯模糊的權重
   static std::vector<double> gaussianFactor1(double sigma) {
     double sigma2 = squ(sigma);
     size_t width = std::max(ceil((double)(sigma * 2)), 0.0);
@@ -38,6 +36,8 @@ private:
     factor[width] = 1.0;
     return factor;
   }
+  
+  //! 回傳gradiance的權重
   static std::vector<double> gradianceFactor1(double sigma) {
     double sigma2 = squ(sigma), ss = sigma * 2;
     size_t width = std::max(ceil(ss), 1.0);
@@ -49,6 +49,8 @@ private:
     factor[width] = 0.0;
     return factor;
   }
+  
+  //! 針對某一方向用某種權重模糊
   Bitmap xyBlur(std::vector<double> const& factor,
                 ssize_t dx, ssize_t dy) const {
     Bitmap ret(*this);
@@ -128,7 +130,7 @@ public:
   void reset(size_t h, size_t w, Pixel const& p) {
     matrix_.reset(h, w, p);
   }
-  
+
   /*!
    * @brief 清除資料, 寬高階規零
    */
@@ -233,6 +235,28 @@ public:
   }
 
   /*!
+   * @brief 回傳矩陣形式
+   */
+  Matrix<Pixel> const& matrix() const {
+    return matrix_;
+  }
+
+  /*!
+   * @brief 回傳矩陣形式 (non-constant form)
+   */
+  Matrix<Pixel>& matrixGet() {
+    return matrix_;
+  }
+
+  /*!
+   * @brief 直接設定
+   */
+  Matrix<Pixel> const& matrix(Matrix<Pixel> const& p) {
+    matrix_.copyFrom(p);
+    return matrix();
+  }
+
+  /*!
    * @brief 回傳高斯模糊
    *
    * @param [in] radiusY 高斯模糊的Y軸方向的sigma
@@ -278,7 +302,7 @@ public:
   Bitmap<Pixel>& gradiancedX(double radiusY, double radiusX) {
     return copyFrom(gradianceX(radiusY, radiusX));
   }
-  
+
   /*!
    * @brief 回傳對y偏微分
    *
@@ -322,48 +346,42 @@ public:
   Pixel const& operator()(size_t y, size_t x, Pixel const& p) const {
     return pixel(y, x, p);
   }
-  
+
   /*! @brief 將資料寫入檔案
    *
    *  @note 未完成, 輸入參數 fg 無用
    */
   bool write(FILE* f, bool bin, unsigned int fg) const {
-    size_t w = width(), h = height();
+    if (fg & 1)
+      return false;
     if (bin) {
-      if (fwrite(&h, sizeof(size_t), 1, f) < 1) return false;
-      if (fwrite(&w, sizeof(size_t), 1, f) < 1) return false;
+      long tmp;
+      if (fwrite(&(tmp = matrix_.cols()), sizeof(tmp), 1, f) < 1) return false;
+      if (fwrite(&(tmp = matrix_.rows()), sizeof(tmp), 1, f) < 1) return false;
     }
     else {
-      if (fprintf(f, "%lu %lu\n", h, w) < 2) return false;
-    }
-    if (fg) {
-      // TODO
-      return false;
+      if (fprintf(f, "%ld %ld\n", (long)matrix_.cols(), (long)matrix_.rows())
+          < 2) return false;
     }
     return true;
-    //return propertyWrite(__f, __bin, __fg);
   }
-  
+
   /*! @brief 將資料讀入
    *
    *  @note 未完成, 輸入參數 fg 無用
    */
   bool read(FILE* f, bool bin, unsigned int fg) {
-    size_t w, h;
-    if (bin) {
-      if (fread(&h, sizeof(size_t), 1, f) < 1) return false;
-      if (fread(&w, sizeof(size_t), 1, f) < 1) return false;
-    }
-    else {
-      if (fscanf(f, "%lu %lu\n", &h, &w) < 2) return false;
-    }
-    if (fg) {
-      // TODO
+    if (fg & 1)
       return false;
+    long tmp1, tmp2;
+    if (bin) {
+      if (fread(&tmp1, sizeof(tmp1), 1, f) < 1) return false;
+      if (fread(&tmp2, sizeof(tmp2), 1, f) < 1) return false;
     }
     else {
-      reset(h, w, Pixel(0));
+      if (fscanf(f, "%ld %ld", &tmp1, &tmp2) < 2) return false;
     }
+    matrix_.size(tmp1, tmp2, Pixel(0));
     return true;
   }
 
@@ -374,7 +392,7 @@ public:
   ObjBase* create() const {
     return new Bitmap();
   }
-  
+
   /*! @brief 複製資料
    *
    * 輸入型別是 \c ObjBase \c const*
@@ -387,16 +405,15 @@ public:
   ObjBase* copyFrom(ObjBase const* b) {
     return &(copyFrom(*(Bitmap*)b));
   }
-  
+
   /*! @brief 回傳class的type
    *
    *  @return \c char \c const\c * 形式的typename
    */
   char const* ctype() const{
-    static char const* ptr = typeid(*this).name();
-    return ptr;
+    return typeid(*this).name();
   }
-  
+
   /*! @brief 回傳class的type
    *
    *  @return \c std::string 形式的typename
@@ -406,6 +423,6 @@ public:
   }
 };
 
-}
+} // meow
 
 #endif // gra_Bitmap_H__
diff --git a/meowpp/gra/BundleAdjustment.h b/meowpp/gra/BundleAdjustment.h
new file mode 100644
index 0000000..9dd01b7
--- /dev/null
+++ b/meowpp/gra/BundleAdjustment.h
@@ -0,0 +1,54 @@
+#ifndef   BundleAdjustment_H__
+#define   BundleAdjustment_H__
+
+#include "Eye.h"
+
+#include "../oo/ObjBase.h"
+
+namespace meow {
+
+enum SceneInfoFlags {
+  CAN_OFFSET = 0x01,
+  CAN_ROTATE = 0x02,
+  CAN_ZOOM   = 0x04
+};
+
+template<class Pixel>
+struct SceneInfo {
+  Eye<Pixel>*    eye;
+  unsigned long flag;
+  
+  SceneInfo(): eye(NULL), flag(0) {
+  }
+  
+  SceneInfo(Eye<Pixel>* e, unsigned long f): eye(e), flag(f) {
+  }
+  
+  SceneInfo(SceneInfo const& si): eye(si.eye), flag(si.flag) {
+  }
+  
+  ~SceneInfo() {
+  }
+};
+
+template<class Pixel>
+class BundleAdjustment: public ObjBase {
+protected:
+  BundleAdjustment() {
+  }
+public:
+  virtual ~BundleAdjustment() {
+  }
+  
+  virtual bool adjustEye(std::vector<SceneInfo<Pixel> >* seq) const {
+    return false;
+  }
+  
+  virtual bool adjustFixedPoint(std::vector<SceneInfo<Pixel> >* seq) const {
+    return false;
+  }
+};
+
+} // meow
+
+#endif // BundleAdjustment_H__
diff --git a/meowpp/gra/BundleAdjustment_LM.h b/meowpp/gra/BundleAdjustment_LM.h
new file mode 100644
index 0000000..7a86666
--- /dev/null
+++ b/meowpp/gra/BundleAdjustment_LM.h
@@ -0,0 +1,375 @@
+#ifndef   BundleAdjustment_LM_H__
+#define   BundleAdjustment_LM_H__
+
+#include "Eye.h"
+
+#include "BundleAdjustment.h"
+
+#include "../math/methods.h"
+#include "../math/Vector.h"
+#include "../math/Matrix.h"
+#include "../math/utility.h"
+
+#include "../oo/ObjBase.h"
+
+#include <algorithm>
+
+namespace meow {
+
+template<class Pixel>
+class BundleAdjustment_LM: public BundleAdjustment<Pixel> {
+private:
+  struct MatchPair {
+    size_t         from_i;
+    Matrix<double> from_m;
+    size_t         to_i;
+    Matrix<double> to_m;
+  };
+  
+  class NoOffsetController {
+  private:
+    std::vector<MatchPair               >* pairs_;
+    std::vector<Rotation3D     <double> >* rot_;
+    std::vector<PhotoProjection<double> >* pho_;
+    BallProjection<double>                 bll_;
+    bool flag;
+  public:
+    NoOffsetController(std::vector<MatchPair               >* pairs,
+                       std::vector<Rotation3D<double>      >* rot,
+                       std::vector<PhotoProjection<double> >* pho,
+                       bool fg):
+    pairs_(pairs), rot_(rot), pho_(pho), bll_(3.0, 1.0), flag(fg) {
+    }
+    Vector<double> init() {
+      if (flag == false) {
+        Vector<double> v(4, 0.0);
+        int i = (*pho_).size() - 1, n = 0;
+        v.entry(n++, (*pho_)[i].focal());
+        v.entry(n++, (*rot_)[i].theta(0));
+        v.entry(n++, (*rot_)[i].theta(1));
+        v.entry(n++, (*rot_)[i].theta(2));
+        return v;
+      }
+      Vector<double> v((*rot_).size() * 4, 0.0);
+      for (size_t n = 0, i = 0, I = (*rot_).size(); i < I; ++i) {
+        v.entry(n++, (*pho_)[i].focal());
+        v.entry(n++, (*rot_)[i].theta(0));
+        v.entry(n++, (*rot_)[i].theta(1));
+        v.entry(n++, (*rot_)[i].theta(2));
+      }
+      return v;
+    }
+    Vector<double> residure(Vector<double> const& v) {
+      for (size_t n = 0, i = 0, I = (*rot_).size(); i < I; ++i) {
+        if (flag == false && i != I - 1) continue;
+        (*pho_)[i].focal(v(n++));
+        (*rot_)[i].theta(0, v(n++));
+        (*rot_)[i].theta(1, v(n++));
+        (*rot_)[i].theta(2, v(n++));
+      }
+      for (typename std::vector<MatchPair>::iterator
+           it = (*pairs_).begin(), ed = (*pairs_).end(); it != ed; ++it) {
+        (*it).from_m.entry(2, 0, -(*pho_)[(*it).from_i].focal());
+        (*it).  to_m.entry(2, 0, -(*pho_)[(*it).  to_i].focal());
+      }
+      //
+      Vector<double> ret((*pairs_).size() * 2, 0.0);
+      size_t n = 0;
+      for (typename std::vector<MatchPair>::iterator
+           it = (*pairs_).begin(), ed = (*pairs_).end(); it != ed; ++it, ++n) {
+        Matrix<double> tr(
+          (*it).to_m
+          -
+          (*pho_)[(*it).to_i].transformate(
+            (*rot_)[(*it).to_i].transformate(
+              (*rot_)[(*it).from_i].transformateInv(
+                bll_.transformate(
+                  (*it).from_m
+                )
+              )
+            )
+          )
+        );
+        ret.entry(n * 2    , tr(0, 0));
+        ret.entry(n * 2 + 1, tr(1, 0));
+      }
+      //*
+      //for (size_t i = 0; i < ret.dimension(); i++) {
+      //  printf("(%-9.1f ", ret(i));
+      //  if (i % 8 == 7) printf("\n");
+      //}
+      //printf("\n");
+      printf("re = %20.5f ?? %f\n", ret.length2(), 25.0 * (*pairs_).size());
+      // */
+      return ret;
+    }
+    Matrix<double> jacobian() const {
+      Matrix<double> ret;
+      if (flag == false) {
+        ret.reset((*pairs_).size() * 2, 4, 0.0);
+      }
+      else {
+        ret.reset((*pairs_).size() * 2, (*rot_).size() * 4, 0.0);
+      }
+      size_t n = 0;
+      for (typename std::vector<MatchPair>::iterator
+           it = (*pairs_).begin(), ed = (*pairs_).end(); it != ed; ++it, ++n) {
+        Matrix<double> tr_A(  bll_               .transformate((*it).from_m));
+        Matrix<double> tr_B((*rot_)[(*it).from_i].transformateInv(tr_A));
+        Matrix<double> tr_C((*rot_)[(*it).  to_i].transformate   (tr_B));
+        Matrix<double> ja_A((*pho_)[(*it).  to_i].jacobian   (tr_C));
+        Matrix<double> ja_B((*rot_)[(*it).  to_i].jacobian   (tr_B));
+        Matrix<double> ja_C((*rot_)[(*it).from_i].jacobianInv(tr_A));
+        Matrix<double> m;
+        int fr = (flag == false ? 0 : (*it).from_i);
+        int to = (flag == false ? 0 : (*it).  to_i);
+        if (flag != false || (*it).from_i == (*rot_).size() - 1) {
+          m = ja_A * ja_B * ja_C * bll_.jacobian((*it).from_m).col(2);
+          ret.entry(n * 2    , fr * 4, m(0, 0));
+          ret.entry(n * 2 + 1, fr * 4, m(1, 0));
+          for (size_t k = 0; k < 3; ++k) {
+            m = ja_A * ja_B * (*rot_)[(*it).from_i].jacobianInv(tr_A, k);
+            ret.entry(n * 2    , fr * 4 + 1 + k, m(0, 0));
+            ret.entry(n * 2 + 1, fr * 4 + 1 + k, m(1, 0));
+          }
+        }
+        if (flag != false || (*it).to_i == (*rot_).size() - 1) {
+          m = (*pho_)[(*it).to_i].jacobian(tr_C, 0);
+          ret.entry(n * 2    , to * 4, m(0, 0));
+          ret.entry(n * 2 + 1, to * 4, m(1, 0));
+          for (size_t k = 0; k < 3; ++k) {
+            m = ja_A * (*rot_)[(*it).to_i].jacobian(tr_B, k);
+            ret.entry(n * 2    , to * 4 + 1 + k, m(0, 0));
+            ret.entry(n * 2 + 1, to * 4 + 1 + k, m(1, 0));
+          }
+        }
+      }
+      return ret;
+    }
+    Matrix<double> identity() const {
+      if (flag == false) {
+        Matrix<double> ret(4, 4, 0.0);
+        ret.identitied();
+        return ret;
+      }
+      Matrix<double> ret((*rot_).size() * 4, (*rot_).size() * 4, 0.0);
+      ret.identitied();
+      return ret;
+    }
+  };
+
+  struct Myself {
+    double t_;
+
+    Myself(): t_(5.0) {
+    }
+
+    Myself(Myself const& m): t_(m.t_) {
+    }
+
+    ~Myself() {
+    }
+  };
+
+  Self<Myself> const self;
+
+public:
+  BundleAdjustment_LM(): self() {
+  }
+
+  BundleAdjustment_LM(BundleAdjustment_LM const& b):
+  self(b.self, Self<Myself>::COPY_FROM) {
+  }
+
+  ~BundleAdjustment_LM() {
+  }
+
+  BundleAdjustment_LM& copyFrom(BundleAdjustment_LM const& b) {
+    self().copyFrom(b.self);
+    return *this;
+  }
+
+  BundleAdjustment_LM& referenceFrom(BundleAdjustment_LM const& b) {
+    self().referenceFrom(b.self);
+    return *this;
+  }
+
+  double threshold() const {
+    return self->t_;
+  }
+
+  double threshold(double t) {
+    self()->t_ = t;
+    return threshold();
+  }
+
+  bool adjustEye(std::vector<SceneInfo<Pixel> >* seq) const {
+    // check
+    size_t N = (*seq).size();
+    for (size_t i = 0; i < N; ++i) {
+      if (  (*seq)[i].flag & CAN_OFFSET ) return false;
+      if (!((*seq)[i].flag & CAN_ROTATE)) return false;
+      if (!((*seq)[i].flag & CAN_ZOOM  )) return false;
+    }
+    // get all pairs
+    typedef typename Camera<Pixel>::FixedPoints2D::IdentityPointsMapIterK FPS_K;
+    std::vector<FPS_K> beg(N), end(N);
+    FPS_K it1, it2;
+    for (size_t i = 0; i < N; ++i) {
+      beg[i] = (*seq)[i].eye->camera().fixedPoints2D().identityPoints().begin();
+      end[i] = (*seq)[i].eye->camera().fixedPoints2D().identityPoints().end  ();
+    }
+    std::vector<std::vector<std::vector<MatchPair> > > all_pairs(N);
+    for (size_t i = 0; i < N; ++i) {
+      all_pairs[i].resize(N);
+      for (size_t j = 0; j < N; ++j) {
+        if (i == j) continue;
+        for (it1 = beg[i], it2 = beg[j]; it1 != end[i] && it2 != end[j]; ) {
+          if     (it1->first < it2->first) ++it1;
+          else if(it1->first > it2->first) ++it2;
+          else {
+            MatchPair tmp;
+            tmp.from_i = i;
+            tmp.from_m = it1->second.matrix();
+            tmp.from_m.rows(3, -(*seq)[i].eye->camera().photo().focal());
+            tmp.  to_i = j;
+            tmp.  to_m = it2->second.matrix();
+            tmp.  to_m.rows(3, -(*seq)[j].eye->camera().photo().focal());
+            all_pairs[i][j].push_back(tmp);
+            ++it1;
+            ++it2;
+          }
+        }
+      }
+    }
+    //
+    std::vector<MatchPair               > pairs;
+    std::vector<Rotation3D     <double> > rot;
+    std::vector<PhotoProjection<double> > pho;
+    std::vector<int                     > sum (N, 0);
+    std::vector<int                     > big (N, 0);
+    std::vector<int                     > real(N);
+    std::vector<int                     > rev (N);
+    // first!!
+    int best = -1, best_ct = -1;
+    for (size_t i = 0; i < N; ++i) {
+      int d = 0;
+      for (size_t j = 0; j < N; ++j) {
+        d += all_pairs[i][j].size();
+      }
+      if (best_ct < d) {
+        best_ct = d;
+        best = i;
+      }
+      real[i] = i;
+      big [i] = i;
+    }
+    for (size_t i = 0; i < N; ++i) {
+      // update else
+      real[i] = best;
+      rev[best] = i;
+      for (size_t j = 0; j < N; ++j) {
+        if ((int)j == best) {
+          sum[j] = -1;
+          continue;
+        }
+        if (sum[j] < 0) continue;
+        sum[j] += all_pairs[best][j].size();
+        if (all_pairs[j][big[j]].size() < all_pairs[j][best].size())
+          big[j] = best;
+      }
+      // add me
+      printf("%d  same as %d\n", best, big[best]);
+      if (big[best] == best) {
+        rot.push_back((*seq)[big[best]].eye->camera()        .rotation());
+        pho.push_back((*seq)[big[best]].eye->camera().photo().projection());
+      }
+      else {
+        rot.push_back(rot[rev[big[best]]]);
+        pho.push_back(pho[rev[big[best]]]);
+      }
+      for (size_t j = 0; j < N; ++j) {
+        if (sum[j] >= 0) continue;
+        for (size_t k = 0, K = all_pairs[best][j].size(); k < K; ++k) {
+          all_pairs[best][j][k].from_i = rev[all_pairs[best][j][k].from_i];
+          all_pairs[best][j][k].  to_i = rev[all_pairs[best][j][k].  to_i];
+          pairs.push_back(all_pairs[best][j][k]);
+        }
+        for (size_t k = 0, K = all_pairs[j][best].size(); k < K; ++k) {
+          all_pairs[j][best][k].from_i = rev[all_pairs[j][best][k].from_i];
+          all_pairs[j][best][k].  to_i = rev[all_pairs[j][best][k].  to_i];
+          pairs.push_back(all_pairs[j][best][k]);
+        }
+      }
+      for (size_t j = 0; j < N; ++j) { printf("%4d ",  sum[j]); } printf("\n");
+      for (size_t j = 0; j < N; ++j) { printf("%4d ",  big[j]); } printf("\n");
+      for (size_t j = 0; j < N; ++j) { printf("%4d ", real[j]); } printf("\n");
+      printf("\n");
+      // bundle adjustment
+      if (i > 0) {
+        NoOffsetController ct1(&pairs, &rot, &pho, false);
+        ct1.residure(levenbergMarquardtTraining(ct1, ct1.init(),
+                                                7.3, 1.0,
+                                                squ(threshold()) * pairs.size(),
+                                                7, 1000));
+        /*
+        NoOffsetController ct2(&pairs, &rot, &pho, true);
+        ct2.residure(levenbergMarquardtTraining(ct2, ct2.init(),
+                                                0.1 * pairs.size(), 1.0,
+                                                squ(threshold()) * pairs.size(),
+                                                7, 1000));
+        // */
+      }
+      // find next
+      best_ct = 0;
+      for (size_t j = 0; j < N; ++j) {
+        if (best_ct < sum[j]) {
+          best    = j;
+          best_ct = sum[j];
+        }
+      }
+      getchar();
+      //if (i == N - 2) break;
+    }
+    // feedback
+    for (size_t i = 0; i < N; ++i) {
+      int id = real[i];
+      (*seq)[id].eye->cameraGet()           .rotation  (rot[i]);
+      (*seq)[id].eye->cameraGet().photoGet().projection(pho[i]);
+      //if (i == N - 2) break;
+    }
+    return true;
+  }
+
+  bool adjustFixedPoint(std::vector<SceneInfo<Pixel> >* seq) const {
+    return false;
+  }
+
+  bool write(FILE* f, bool bin, unsigned int fg) const {
+    return false;
+  }
+
+  bool read(FILE* f, bool bin, unsigned int fg) const {
+    return false;
+  }
+
+  ObjBase* create() const {
+    return new BundleAdjustment_LM;
+  }
+
+  ObjBase* copyFrom(ObjBase const* o) {
+    return &(copyFrom(*(BundleAdjustment_LM const*)o));
+  }
+
+  char const* ctype() const {
+    return typeid(*this).name();
+  }
+
+  std::string type() const {
+    return std::string(ctype());
+  }
+};
+
+} // meow
+
+#endif // BundleAdjustment_LM_H__
diff --git a/meowpp/gra/Camera.h b/meowpp/gra/Camera.h
index 45321de..2b6347c 100644
--- a/meowpp/gra/Camera.h
+++ b/meowpp/gra/Camera.h
@@ -32,13 +32,10 @@ private:
     Myself() {
       fixed2D_.dimension(2);
     }
-    ~Myself() {
+    Myself(Myself const& v):
+    photo_(v.photo_), rot_(v.rot_), fixed2D_(v.fixed2D_) {
     }
-    Myself& copyFrom(Myself const& b) {
-      photo_  .copyFrom(b.  photo_);
-      rot_    .copyFrom(b.    rot_);
-      fixed2D_.copyFrom(b.fixed2D_);
-      return *this;
+    ~Myself() {
     }
   };
 
@@ -47,14 +44,13 @@ public:
   /*!
    * @brief constructor
    */
-  Camera(): self(true) {
+  Camera(): self() {
   }
   
   /*!
    * @brief copy constructor
    */
-  Camera(Camera const& b): self(false) {
-    copyFrom(b);
+  Camera(Camera const& b): self(b.self, Self<Myself>::COPY_FROM) {
   }
 
   /*!
@@ -119,7 +115,7 @@ public:
    * @brief 設定rotation
    */
   Rotation3D<double> const& rotation(Rotation3D<double> const& rot) {
-    self()->rot_ = rot;
+    self()->rot_.copyFrom(rot);
     return rotation();
   }
 
@@ -142,7 +138,7 @@ public:
    */
   FixedPoints2D const& fixedPoints2D(FixedPoints2D const& fps2d) const {
     if (fps2d.dimension() == 2) {
-      self()->fixed2D_ = fps2d;
+      self()->fixed2D_.copyFrom(fps2d);
     }
     return fixedPoints2D();
   }
@@ -150,7 +146,7 @@ public:
   /*!
    * @brief 取得編號為i的fixed points 2d
    */
-  Vector<double> fixedPoints2D(int i) {
+  Vector<double> fixedPoint2D(int i) {
     return self->fixed2D_.identityPoint(i);
   }
 
@@ -159,7 +155,7 @@ public:
    */
   bool inside(Vector3D<double> p) const {
     return self->photo_.inside(
-      Vector3D<double>(self->rot_.transformate(p.matrix())));
+      Vector3D<double>(rotation().transformate(p.matrix())));
   }
 
   /*!
@@ -167,7 +163,7 @@ public:
    */
   Pixel color(Vector3D<double> p) const {
     return self->photo_.color(
-      Vector3D<double>(self->rot_.transformate(p.matrix())));
+      Vector3D<double>(rotation().transformate(p.matrix())));
   }
   
   /*!
@@ -182,7 +178,20 @@ public:
    *  @note 未完成
    */
   bool write(FILE* f, bool bin, unsigned int fg) const {
-    return false;
+    if (bin) {
+      double tmp;
+      for (size_t i = 0; i < 3; ++i) {
+        if (fwrite(&(tmp = rotation().theta(i)), sizeof(tmp), 1, f) < 1)
+          return false;
+      }
+    }
+    else {
+      for (size_t i = 0; i < 3; ++i) {
+        if (fprintf(f, "%f ", rotation().theta(i)) < 1) return false;
+      }
+      fprintf(f, "\n");
+    }
+    return (fixedPoints2D().write(f, bin, fg) && photo().write(f, bin, fg));
   }
 
   /*! @brief 將資料讀入
@@ -190,7 +199,23 @@ public:
    *  @note 未完成
    */
   bool read(FILE* f, bool bin, unsigned int fg) {
-    return false;
+    if (bin) {
+      double tmp;
+      for (size_t i = 0; i < 3; ++i) {
+        if (fread(&tmp, sizeof(tmp), 1, f) < 1) {
+          return false;
+        }
+        rotationGet().theta(i, tmp);
+      }
+    }
+    else {
+      double a;
+      for (size_t i = 0; i < 3; ++i) {
+        if (fscanf(f, "%lf", &a) < 1) return false;
+        rotationGet().theta(i, a);
+      }
+    }
+    return (fixedPoints2DGet().read(f, bin, fg) && photoGet().read(f, bin, fg));
   }
 
   /*! @brief new一個自己
@@ -229,316 +254,8 @@ public:
   std::string type() const {
     return std::string(ctype());
   }
-  
-  ////////////////////////////////////////////////////////////////////
-private:
-  class BoundleAdjustment2D {
-  private:
-    class Parameters {
-    private:
-      std::vector<Camera>&                  cam_;
-      std::vector<Rotation3D<double> >      rot_;
-      std::vector<PhotoProjection<double> > pho_;
-      struct Pair {
-        size_t i1_;
-        size_t i2_;
-        Vector<double> v1_;
-        Vector<double> v2_;
-        Pair(size_t a, size_t b,
-             Vector<double> const& v1, Vector<double> const& v2):
-        i1_(a), i2_(b), v1_(v1), v2_(v2) {
-        }
-      };
-      std::vector<Pair> pairs_;
-
-      void setParameters(Vector<double> const& v) {
-        size_t n = 0;
-        for (size_t i = 0, I = cam_.size(); i < I; ++i) {
-          pho_[i].focal(v(n++));
-          for (size_t j = 0; j < 3; j++) {
-            rot_[i].theta(j, v(n++));
-          }
-        }
-        for (size_t i = 0, I = pairs_.size(); i < I; ++i) {
-          pairs_[i].v1_.entry(2, pho_[pairs_[i].i1_].focal());
-          pairs_[i].v2_.entry(2, pho_[pairs_[i].i2_].focal());
-        }
-      }
-      Vector<double> getParameters() const {
-        Vector<double> ret(cam_.size() * 4, 0.0);
-        for (size_t i = 0, I = cam_.size(); i < I; ++i) {
-          ret.entry(i * 4, pho_[i].focal());
-          for (size_t j = 0; j < 3; ++j) {
-            ret.entry(i * 4 + 1 + j, rot_[i].theta(j));
-          }
-        }
-        return ret;
-      }
-      Vector<double> residureV() const {
-        Vector<double> ret(pairs_.size() * 3, 0.0);
-        for (size_t i = 0, I = pairs_.size(); i < I; ++i) {
-          size_t i_from = pairs_[i].i1_;
-          size_t i_to   = pairs_[i].i2_;
-          Matrix<double> v_from(pairs_[i].v1_.matrix());
-          Matrix<double> v_to  (pairs_[i].v2_.matrix());
-          Matrix<double> v_tr(
-            pho_[i_to].transformate(
-              rot_[i_to].transformate(
-                rot_[i_from].transformateInv(
-                  BallProjection<double>(3, 1.0).transformate(
-                    v_from
-                  )
-                )
-              )
-            )
-          );
-          Matrix<double> delta(v_to - v_tr);
-          for (size_t j = 0; j < 3; ++j) {
-            ret.entry(i * 3 + j, delta(j, 0));
-          }
-        }
-        return ret;
-      }
-    public:
-      Parameters(std::vector<Camera>& cam): cam_(cam) {
-        rot_.resize(cam_.size());
-        pho_.resize(cam_.size(), PhotoProjection<double>(3));
-        for (size_t i = 0, I = cam_.size(); i < I; ++i) {
-          rot_[i].referenceFrom(cam_[i].rotation());
-          pho_[i].focal(cam_[i].photo().focal());
-        }
-        for (size_t i = 0, I = cam_.size(); i < I; ++i) {
-          std::map<int,Vector<double> >const& p1 = (
-            cam_[i].fixedPoints2D().identityPoints());
-          for (size_t j = 0; j < I; ++j) {
-            if (i == j) continue;
-            std::map<int,Vector<double> >const& p2 = (
-              cam_[j].fixedPoints2D().identityPoints());
-            for (std::map<int,Vector<double> >::const_iterator
-                 it1 = p1.begin(); it1 != p1.end(); ++it1) {
-              for (std::map<int,Vector<double> >::const_iterator
-                   it2 = p2.begin(); it2 != p2.end(); ++it2) {
-                if (it1->first != it2->first) continue;
-                Vector<double> v1(it1->second), v2(it2->second);
-                v1.dimension(3, 0.0);
-                v2.dimension(3, 0.0);
-                pairs_.push_back(Pair(i, j, v1, v2));
-              }
-            }
-          }
-        }
-      }
-      Vector<double> init() const {
-        return getParameters();
-      }
-      Vector<double> residure(Vector<double> const& v) const {
-        ((Parameters*)this)->setParameters(v);
-        return residureV();
-      }
-      Matrix<double> jacobian(Vector<double> const& v) const {
-        //setParameters(v);
-        Matrix<double> ret(pairs_.size() * 3, v.dimension(), 0.0);
-        for (size_t i = 0, I = pairs_.size(); i < I; ++i) {
-          for (size_t j = 0, J = v.dimension(); j < J; ++j) {
-            size_t j0 = j / 4;
-            size_t dj = j % 4;
-            size_t i_from = pairs_[i].i1_;
-            size_t i_to   = pairs_[i].i2_;
-            Matrix<double> v_from(pairs_[i].v1_.matrix());
-            Matrix<double> v_to  (pairs_[i].v2_.matrix());
-            Matrix<double> v_tr  (3, 1, 0.0);
-            if (j0 == i_from) {
-              if (dj == 0) {
-                v_tr = (
-                  pho_[i_to].jacobian(
-                    rot_[i_to].transformate(
-                      rot_[i_from].transformateInv(
-                        BallProjection<double>(3, 1.0).transformate(
-                          v_from
-                        )
-                      )
-                    )
-                  )
-                  *
-                  rot_[i_to].jacobian(
-                    rot_[i_from].transformateInv(
-                      BallProjection<double>(3, 1.0).transformate(
-                        v_from
-                      )
-                    )
-                  )
-                  *
-                  rot_[i_from].jacobianInv(
-                    BallProjection<double>(3, 1.0).transformate(
-                      v_from
-                    )
-                  )
-                  *
-                  BallProjection<double>(3, 1.0).jacobian(
-                    v_from
-                  ).col(2)
-                  );
-              }
-              else {
-                v_tr = (
-                  pho_[i_to].jacobian(
-                    rot_[i_to].transformate(
-                      rot_[i_from].transformateInv(
-                        BallProjection<double>(3, 1.0).transformate(
-                          v_from
-                        )
-                      )
-                    )
-                  )
-                  *
-                  rot_[i_to].jacobian(
-                    rot_[i_from].transformateInv(
-                      BallProjection<double>(3, 1.0).transformate(
-                        v_from
-                      )
-                    )
-                  )
-                  *
-                  rot_[i_from].jacobianInv(
-                    BallProjection<double>(3, 1.0).transformate(
-                      v_from
-                    )
-                  )
-                  );
-              }
-            }
-            else if (j0 == i_to) {
-              if (dj == 0) {
-                v_tr = (
-                  pho_[i_to].jacobian(
-                    rot_[i_to].transformate(
-                      rot_[i_from].transformateInv(
-                        BallProjection<double>(3, 1.0).transformate(
-                          v_from
-                        )
-                      )
-                    )
-                  ).col(2)
-                );
-              }
-              else {
-                v_tr = (
-                  pho_[i_to].jacobian(
-                    rot_[i_to].transformate(
-                      rot_[i_from].transformateInv(
-                        BallProjection<double>(3, 1.0).transformate(
-                          v_from
-                        )
-                      )
-                    )
-                  )
-                  *
-                  rot_[i_to].jacobian(
-                    rot_[i_from].transformateInv(
-                      BallProjection<double>(3, 1.0).transformate(
-                        v_from
-                      )
-                    ),
-                    dj - 1
-                  )
-                );
-              }
-            }
-            for (size_t k = 0; k < 3; ++k) {
-              ret.entry(i * 3 + k, j, -v_tr(k, 0));
-            }
-          }
-        }
-        return ret;
-      }
-      Matrix<double> identity(Vector<double> const& v) const {
-        //setParameters(v);
-        Matrix<double> ret(v.dimension(), v.dimension(), 0.0);
-        ret.identity();
-        return ret;
-      }
-      double averageResidure() const {
-        Vector<double> res(residureV());
-        double sum = 0;
-        for (size_t i = 0, I = res.dimension(); i < I; ++i) {
-          sum += res(i);
-        }
-        return sum / res.dimension();
-      }
-      size_t dimensinonI() const {
-        return cam_.size() * 4;
-      }
-      size_t dimensionO() const {
-        return pairs_.size() * 3;
-      }
-    };
-    class F {
-    private:
-      Parameters& p_;
-    public:
-      F(Parameters& p): p_(p) {
-      }
-      Vector<double> operator()(Vector<double> const& v) const {
-        return p_.residure(v);
-      }
-    };
-    class J {
-    private:
-      Parameters& p_;
-    public:
-      J(Parameters& p): p_(p) {
-      }
-      Matrix<double> operator()(Vector<double> const& v) const {
-        return p_.jacobian(v);
-      }
-    };
-    class I {
-    private:
-      Parameters& p_;
-    public:
-      I(Parameters& p): p_(p) {
-      }
-      Matrix<double> operator()(Vector<double> const& v) const {
-        return p_.identity(v);
-      }
-    };
-    class Stop {
-    private:
-      Parameters& p_;
-      double t_;
-    public:
-      Stop(Parameters& p, double t): p_(p), t_(t) {
-      }
-      bool operator()(double r) const {
-        return (r < p_.dimensionO() * t_);
-      }
-    };
-  public:
-    BoundleAdjustment2D() {
-    }
-    ~BoundleAdjustment2D() {
-    }
-    double operator()(std::vector<Camera>* cs, double threshold) const {
-      Parameters p(*cs);
-      Vector<double> v0(p.init());
-      levenbergMarquardt(F(p), J(p), I(p), v0, Stop(p, threshold), 100000);
-      return p.averageResidure();
-    }
-  };
-public:
-  /*!
-   * @brief 將數台camera用fixed points做boundle adjustment
-   *
-   * @param [in] cs        要調整的cameras
-   * @param [in] threshold 允許誤差值
-   * @return 誤差值
-   */
-  static double boundleAdjustment2D(std::vector<Camera>* cs, double threshold) {
-    static BoundleAdjustment2D bdl;
-    return bdl(cs, threshold);
-  }
 };
 
-}
+} // meow
 
 #endif // gra_Camera_H__
diff --git a/meowpp/gra/Eye.h b/meowpp/gra/Eye.h
new file mode 100644
index 0000000..2621a16
--- /dev/null
+++ b/meowpp/gra/Eye.h
@@ -0,0 +1,174 @@
+#ifndef   Eye_H__
+#define   Eye_H__
+
+#include "Camera.h"
+
+#include "../Self.h"
+#include "../oo/ObjBase.h"
+
+namespace meow {
+
+/*!
+ * @brief 一個 \c Camera 加上一個offset transformation
+ * 
+ * @author cat_leopard
+ */
+template<class Pixel>
+class Eye: public ObjBase {
+private:
+  struct Myself {
+    Camera<Pixel>    cam_;
+    Vector3D<double> ofs_;
+    
+    Myself(): cam_(), ofs_(0.0, 0.0, 0.0) {
+    }
+    
+    Myself(Camera<Pixel> const& c, Vector3D<double> const& o): cam_(c), ofs_(o){
+    }
+    
+    Myself(Myself const& b): cam_(b.cam_), ofs_(b.ofs_) {
+    }
+    
+    ~Myself() {
+    }
+  };
+  
+  Self<Myself> const self;
+public:
+  Eye(): self() {
+  }
+  
+  Eye(Eye const& b): self(b.self(), Self<Myself>::COPY_FROM) {
+  }
+  
+  Eye(Camera<Pixel> const& c, Vector3D<double> const& o): self(Myself(c, o)) {
+  }
+  
+  ~Eye() {
+  }
+  
+  Eye& copyFrom(Eye const& e) {
+    self().copyFrom(e.self);
+    return *this;
+  }
+  
+  Eye& referenceFrom(Eye const& e) {
+    self().referenceFrom(e.self);
+    return *this;
+  }
+  
+  Camera<Pixel> const& camera() const {
+    return self->cam_;
+  }
+  
+  Camera<Pixel>& cameraGet() {
+    return self()->cam_;
+  }
+  
+  Camera<Pixel> const& camera(Camera<Pixel> const& c) {
+    self()->cam_.copyFrom(c);
+    return camera();
+  }
+  
+  Vector3D<double> const& offset() const {
+    return self->ofs_;
+  }
+  
+  Vector3D<double>& offsetGet() {
+    return self()->ofs_;
+  }
+  
+  Vector3D<double> const& offset(Vector3D<double> const& ofs) {
+    self()->ofs_ = ofs;
+    return offset();
+  }
+  
+  bool inside(Vector3D<double> const& v) const {
+    return camera().inside(v - offset());
+  }
+  
+  Eye& operator=(Eye const& e) {
+    return copyFrom(e);
+  }
+  
+  /*! @brief 將資料寫入檔案
+   *
+   *  @note 未完成
+   */
+  bool write(FILE* f, bool bin, unsigned int fg) const {
+    if (bin) {
+      double tmp;
+      for (size_t i = 0; i < 3; ++i) {
+        if (fwrite(&(tmp = offset()(i)), sizeof(tmp), 1, f) < 1)
+          return false;
+      }
+    }
+    else {
+      for (size_t i = 0; i < 3; ++i) {
+        if (fprintf(f, "%f ", offset()(i)) < 1) return false;
+      }
+      fprintf(f, "\n");
+    }
+    return camera().write(f, bin, fg);
+  }
+
+  /*! @brief 將資料讀入
+   *
+   *  @note 未完成
+   */
+  bool read(FILE* f, bool bin, unsigned int fg) {
+    if (bin) {
+      double tmp[3];
+      if (fread(tmp, sizeof(double), 3, f) < 3) return false;
+      offsetGet().xyz(tmp[0], tmp[1], tmp[2]);
+    }
+    else {
+      double a, b, c;
+      if (fscanf(f, "%lf %lf %lf", &a, &b, &c) < 3) return false;
+      offsetGet().x(a);
+      offsetGet().y(b);
+      offsetGet().z(c);
+    }
+    return cameraGet().read(f, bin, fg);
+  }
+
+  /*! @brief new一個自己
+   *
+   *  @return 一個new出來的pointer
+   */
+  ObjBase* create() const {
+    return new Eye();
+  }
+
+  /*! @brief 複製資料
+   *
+   * 輸入型別是 \c ObjBase \c const*
+   * 事實上這個method就只是幫忙轉型然後呼叫原本的\c copyFrom
+   *
+   *  @param [in] b 資料來源
+   *  @return this
+   */
+  ObjBase* copyFrom(ObjBase const* b) {
+    return &(copyFrom(*(Eye*)b));
+  }
+
+  /*! @brief 回傳class的type
+   *
+   *  @return \c char \c const\c * 形式的typename
+   */
+  char const* ctype() const{
+    return typeid(*this).name();
+  }
+
+  /*! @brief 回傳class的type
+   *
+   *  @return \c std::string 形式的typename
+   */
+  std::string type() const {
+    return std::string(ctype());
+  }
+};
+
+} // meow
+
+#endif // Eye_H__
diff --git a/meowpp/gra/FeaturePoint.h b/meowpp/gra/FeaturePoint.h
index 4b80a3d..c3ab8f3 100644
--- a/meowpp/gra/FeaturePoint.h
+++ b/meowpp/gra/FeaturePoint.h
@@ -75,6 +75,13 @@ public:
   }
 
   /*!
+   * @brief 回傳position (non-const reference)
+   */
+  Vector<Scalar>& positionGet() {
+    return pos_;
+  }
+
+  /*!
    * @brief 回傳description
    */
   Vector<Description> const& description() const {
@@ -82,6 +89,13 @@ public:
   }
 
   /*!
+   * @brief 回傳description (non-const reference)
+   */
+  Vector<Description>& descriptionGet() {
+    return des_;
+  }
+
+  /*!
    * @brief 修改position
    */
   Vector<Scalar> const& position(Vector<Scalar> const& p) const {
@@ -128,20 +142,6 @@ public:
   }
 
   /*!
-   * @brief 取得position
-   */
-  Vector<Scalar>& positionGet() {
-    return pos_;
-  }
-
-  /*!
-   * @brief 取得description
-   */
-  Vector<Description>& descriptionGet() {
-    return des_;
-  }
-
-  /*!
    * @brief same as copyFrom(fp)
    */
   FeaturePoint& operator=(FeaturePoint const& fp) {
@@ -163,11 +163,53 @@ public:
   }
 
   bool write(FILE* f, bool bin, unsigned int fg) const {
-    return false;
+    if (bin) {
+      double tmp;
+      for (size_t i = 0, I = position().dimension(); i < I; ++i) {
+        if (fwrite(&(tmp = position(i)), sizeof(tmp), 1, f) < 1) return false;
+      }
+      for (size_t i = 0, I = description().dimension(); i < I; ++i) {
+        if (fwrite(&(tmp = description(i)), sizeof(tmp), 1, f) < 1)
+          return false;
+      }
+    }
+    else {
+      for (size_t i = 0, I = position().dimension(); i < I; ++i) {
+        if (fprintf(f, "%f ", (double)position(i)) < 1) return false;
+      }
+      fprintf(f, "\n");
+      for (size_t i = 0, I = description().dimension(); i < I; ++i) {
+        if (fprintf(f, "%f ", (double)description(i)) < 1) return false;
+      }
+      fprintf(f, "\n");
+    }
+    return true;
   }
 
   bool read (FILE* f, bool bin, unsigned int fg) {
-    return false;
+    if (bin) {
+      double tmp;
+      for (size_t i = 0, I = position().dimension(); i < I; ++i) {
+        if (fread(&tmp, sizeof(tmp), 1, f) < 1) return false;
+        position(i, tmp);
+      }
+      for (size_t i = 0, I = description().dimension(); i < I; ++i) {
+        if (fread(&tmp, sizeof(tmp), 1, f) < 1) return false;
+        description(i, tmp);
+      }
+    }
+    else {
+      double tmp;
+      for (size_t i = 0, I = position().dimension(); i < I; ++i) {
+        if (fscanf(f, "%lf", &tmp) < 1) return false;
+        position(i, tmp);
+      }
+      for (size_t i = 0, I = description().dimension(); i < I; ++i) {
+        if (fscanf(f, "%lf", &tmp) < 1) return false;
+        description(i, tmp);
+      }
+    }
+    return true;
   }
 
   ObjBase* create() const {
@@ -179,8 +221,7 @@ public:
   }
 
   char const* ctype() const {
-    static char const* ptr = typeid(*this).name();
-    return ptr;
+    return typeid(*this).name();
   }
 
   std::string type() const {
@@ -188,6 +229,6 @@ public:
   }
 };
 
-}
+} // meow
 
 #endif // gra_FeaturePoint_H__
diff --git a/meowpp/gra/FeaturePointsDetector.h b/meowpp/gra/FeaturePointsDetector.h
index 415419a..858c66a 100644
--- a/meowpp/gra/FeaturePointsDetector.h
+++ b/meowpp/gra/FeaturePointsDetector.h
@@ -18,9 +18,9 @@ public:
   virtual ~FeaturePointsDetector() { }
 
   virtual std::vector<FeaturePoint<double, double> >
-  detect(Bitmap<Pixel> const& __bitmap) const = 0;
+  detect(Bitmap<Pixel> const& bitmap) const = 0;
 };
 
-}
+} // meow
 
 #endif // gra_FeaturePointsDetector_H__
diff --git a/meowpp/gra/FeaturePointsDetector_Harris.h b/meowpp/gra/FeaturePointsDetector_Harris.h
index 25a46b2..11529f0 100644
--- a/meowpp/gra/FeaturePointsDetector_Harris.h
+++ b/meowpp/gra/FeaturePointsDetector_Harris.h
@@ -13,7 +13,6 @@
 
 #include <vector>
 
-
 namespace meow {
 
 /*!
@@ -33,27 +32,26 @@ private:
     double     lightL_;
     double   featureG_;
     size_t     boundB_;
-    
-    Myself() {
-      ratioK_     = 0.03;
-      thresholdR_ = 0.001;
-      sizeW_      = 2.0;
-      noiseN_     = 3.0;
-      lightL_     = 30.0;
-      featureG_   = 3.0;
-      boundB_     = 10u;
+
+    Myself():
+    ratioK_(0.03),
+    thresholdR_(0.001),
+    sizeW_(2.0),
+    noiseN_(3.0),
+    lightL_(30.0),
+    featureG_(3.0),
+    boundB_(10u) {
     }
-    ~Myself() {
+    Myself(Myself const& m):
+    ratioK_(m.ratioK_),
+    thresholdR_(m.thresholdR_),
+    sizeW_(m.sizeW_),
+    noiseN_(m.noiseN_),
+    lightL_(m.lightL_),
+    featureG_(m.featureG_),
+    boundB_(m.boundB_){
     }
-    Myself& copyFrom(Myself const& b) {
-      ratioK_     = b.ratioK_     ;
-      thresholdR_ = b.thresholdR_ ;
-      sizeW_      = b.sizeW_      ;
-      noiseN_     = b.noiseN_     ;
-      lightL_     = b.lightL_     ;
-      featureG_   = b.featureG_   ;
-      boundB_     = b.boundB_     ;
-      return *this;
+    ~Myself() {
     }
   };
 
@@ -62,19 +60,11 @@ public:
   typedef FeaturePoint<double, double> MyFeaturePoint;
   typedef std::vector<MyFeaturePoint>  MyFeaturePoints;
   //! @brief constructor 使用預設參數
-  FPD_Harris(): self(true) {
-    self()->ratioK_     = 0.03;
-    self()->thresholdR_ = 0.001;
-    self()->sizeW_      = 2.0;
-    self()->noiseN_     = 3.0;
-    self()->lightL_     = 30.0;
-    self()->featureG_   = 3.0;
-    self()->boundB_     = 10u;
+  FPD_Harris(): self() {
   }
   
   //! @brief constructor 參數複製自另一個 FeaturePointsDetector_Harris
-  FPD_Harris(FPD_Harris const& fps): self(false) {
-    self().copyFrom(fps.self);
+  FPD_Harris(FPD_Harris const& fps): self(fps.self, Self<Myself>::COPY_FROM) {
   }
   
   //! @brief 解構子
@@ -178,9 +168,11 @@ public:
   MyFeaturePoints detect(Bitmap<Pixel> const& bmp) const {
     Bitmap<Pixel> input = bmp;
     
+    // gradiance
     Bitmap<Pixel> input_gx(input.gradianceX(0, self->noiseN_));
     Bitmap<Pixel> input_gy(input.gradianceY(self->noiseN_, 0));
     
+    // get Matrix I for each pixel
     Bitmap<double> Ixx(input.height(), input.width(), 0.0);
     Bitmap<double> Iyy(input.height(), input.width(), 0.0);
     Bitmap<double> Ixy(input.height(), input.width(), 0.0);
@@ -194,10 +186,12 @@ public:
       }
     }
     
+    // blur
     Ixx.gaussianed(self->sizeW_, self->sizeW_);
     Iyy.gaussianed(self->sizeW_, self->sizeW_);
     Ixy.gaussianed(self->sizeW_, self->sizeW_);
     
+    // filter too flat or on edge
     Bitmap<double>  R(input.height(), input.width(),   0.0);
     Bitmap<bool> good(input.height(), input.width(), false);
     ssize_t b = self->boundB_;
@@ -211,6 +205,7 @@ public:
       }
     }
     
+    // find union neighbor
     DisjointSet dsj(input.size());
     ssize_t dy[2] = {0, 1};
     ssize_t dx[2] = {1, 0};
@@ -227,6 +222,7 @@ public:
       }
     }
 
+    // find local maximum
     std::vector<size_t> max_i(input.size());
     for (size_t i = 0, I = input.size(); i < I; i++) {
       max_i[i] = i;
@@ -239,6 +235,7 @@ public:
       }
     }
     
+    // blur before get description
     input.gaussianed(self->featureG_, self->featureG_);
 
     MyFeaturePoints ret;
@@ -253,7 +250,7 @@ public:
         }
         ssize_t dx[4] = {1, 0, -1, 0};
         ssize_t dy[4] = {0, 1, 0, -1};
-        std::vector<double> desc;
+        std::vector<double> desc; // description
         for (ssize_t d = 1; d <= (ssize_t)self->boundB_; d++) {
           std::vector<double> light;
           size_t max_id = 0;
diff --git a/meowpp/gra/FeaturePointsMatch.h b/meowpp/gra/FeaturePointsMatch.h
index 8b05632..ffdbe53 100644
--- a/meowpp/gra/FeaturePointsMatch.h
+++ b/meowpp/gra/FeaturePointsMatch.h
@@ -3,32 +3,15 @@
 
 #include "FeaturePoint.h"
 
+#include "../utility.h"
 #include "../oo/ObjBase.h"
 
-#include <utility>
 #include <cstdlib>
 
 namespace meow {
 
-struct FeaturePointIndexPair {
-  std::pair<size_t, size_t> from;
-  std::pair<size_t, size_t> to;
-  
-  FeaturePointIndexPair() {
-  }
-  FeaturePointIndexPair(size_t ff, size_t fs,
-                        size_t tf, size_t ts) {
-    from.first  = ff;
-    from.second = fs;
-    to.first    = tf;
-    to.second   = ts;
-  }
-  bool operator==(FeaturePointIndexPair const& p) const {
-    return (from == p.from && to == p.to);
-  }
-};
-
-typedef std::vector<FeaturePointIndexPair> FeaturePointIndexPairs;
+typedef PairToPair<size_t, size_t, size_t, size_t> FeaturePointIndexPair ;
+typedef std::vector<FeaturePointIndexPair>         FeaturePointIndexPairs;
 
 template<class Scalar, class Description>
 class FeaturePointsMatch: public ObjBase {
@@ -58,6 +41,6 @@ public:
                                        FeaturePointss const& fpss) const = 0;
 };
 
-}
+} // meow
 
 #endif // gra_FeaturePointsMatch_H__
diff --git a/meowpp/gra/FeaturePointsMatch_K_Match.h b/meowpp/gra/FeaturePointsMatch_K_Match.h
index ca47d27..f1fff12 100644
--- a/meowpp/gra/FeaturePointsMatch_K_Match.h
+++ b/meowpp/gra/FeaturePointsMatch_K_Match.h
@@ -1,13 +1,10 @@
 #ifndef   gra_FeaturePointsMatch_K_Match_H__
 #define   gra_FeaturePointsMatch_K_Match_H__
 
-#include "../dsa/VP_Tree.h"
-
 #include "FeaturePointsMatch.h"
 
 #include "../Self.h"
-
-
+#include "../dsa/VP_Tree.h"
 #include "../oo/ObjBase.h"
 
 #include <cstdlib>
@@ -19,20 +16,20 @@ class FeaturePointsMatch_K_Match:
 public FeaturePointsMatch<Scalar, Description> {
 # define FPMKM FeaturePointsMatch_K_Match
 public:
-  typedef std::vector<FeaturePoint<Scalar, Description> > FeaturePoints;
+  typedef std::vector<FeaturePoint<Scalar, Description> > FeaturePoints ;
   typedef std::vector<FeaturePoints                     > FeaturePointss;
 private:
   struct Node {
-    size_t id_;
-    size_t index_;
-    FeaturePointss const* ptr_;
+    size_t                   id_;
+    size_t                index_;
+    FeaturePointss const*   ptr_;
 
     Node() {
     }
     Node(Node const& nd) {
-      id_    = nd.id_;
+      id_    = nd.   id_;
       index_ = nd.index_;
-      ptr_   = nd.ptr_;
+      ptr_   = nd.  ptr_;
     }
     Node(size_t id, size_t index, FeaturePointss const* ptr) {
       id_    = id;
@@ -48,30 +45,30 @@ private:
       return (*ptr_)[id_][index_][id];
     }
   };
+  
   struct Myself {
     size_t k_;
     Myself() {
       k_ = 1;
     }
-    ~Myself() {
+    Myself(size_t k): k_(k) {
     }
-    Myself& copyFrom(Myself const& m) {
-      k_ = m.k_;
-      return *this;
+    Myself(Myself const& m): k_(m.k_) {
+    }
+    ~Myself() {
     }
   };
 
   Self<Myself> const self;
 public:
-  FPMKM(): self(true) {
+  FPMKM(): self() {
   }
 
-  FPMKM(FPMKM const& m): self(false) {
+  FPMKM(FPMKM const& m): self(m.self, Self<Myself>::COPY_FROM) {
     self().copyFrom(m.self);
   }
 
-  FPMKM(size_t k): self(true) {
-    self()->k_ = k;
+  FPMKM(size_t k): self(Myself(k)) {
   }
 
   ~FPMKM() {
@@ -175,8 +172,7 @@ public:
   }
 
   char const* ctype() const {
-    static char const* ptr = typeid(*this).name();
-    return ptr;
+    return typeid(*this).name();
   }
 
   std::string type() const {
@@ -185,6 +181,6 @@ public:
 # undef FPMKM
 };
 
-}
+} // meow
 
 #endif // gra_FeaturePointsMatch_K_Match_H__
diff --git a/meowpp/gra/IdentityPoints.h b/meowpp/gra/IdentityPoints.h
index c33eb55..1646bf7 100644
--- a/meowpp/gra/IdentityPoints.h
+++ b/meowpp/gra/IdentityPoints.h
@@ -2,12 +2,11 @@
 #define   gra_IdentityPoints_H__
 
 #include "../Self.h"
-
 #include "../math/Vector.h"
-
 #include "../oo/ObjBase.h"
 
 #include <map>
+#include <set>
 
 #include <cstdlib>
 
@@ -24,45 +23,42 @@ public:
   typedef typename std::map<ID, Vector<Scalar> >     IdentityPointsMap;
   typedef typename IdentityPointsMap::      iterator IdentityPointsMapIter;
   typedef typename IdentityPointsMap::const_iterator IdentityPointsMapIterK;
-  
+
 private:
   struct Myself {
     IdentityPointsMap    points_;
     size_t            dimension_;
-    
+
     Myself() {
       dimension_ = 1;
     }
-    ~Myself() {
+    Myself(Myself const& m): points_(m.points_), dimension_(m.dimension_) {
     }
-    Myself& copyFrom(Myself const& b) {
-      points_    = b.points_;
-      dimension_ = b.dimension_;
-      return *this;
+    ~Myself() {
     }
   };
-  
+
   Self<Myself> const self;
 public:
   /*!
    * @brief constructor
    */
-  IdentityPoints(): self(true) {
+  IdentityPoints(): self() {
   }
-  
+
   /*!
    * @brief constructor, 並且複製資料
    */
-  IdentityPoints(IdentityPoints const& b): self(false) {
-    copyFrom(b);
+  IdentityPoints(IdentityPoints const& b):
+  self(b.self, Self<Myself>::COPY_FROM) {
   }
-  
+
   /*!
    * @brief destructor
    */
   ~IdentityPoints() {
   }
-  
+
   /*!
    * @brief 複製資料
    */
@@ -70,7 +66,7 @@ public:
     self().copyFrom(b.self);
     return *this;
   }
-  
+
   /*!
    * @brief 參照
    */
@@ -78,42 +74,42 @@ public:
     self().referenceFrom(b.self);
     return *this;
   }
-  
+
   /*!
    * @brief 清除一切identity points
    */
   void clear() {
     self()->points_.clear();
   }
-  
+
   /*!
    * @brief 回傳有幾個identity points
    */
   size_t size() const {
     return self->points_.size();
   }
-  
+
   /*!
    * @brief 回傳是否沒有identity points
    */
   bool empty() const {
     return (size() == 0u);
   }
-  
+
   /*!
    * @brief 檢查某id是否有使用
    */
   bool exist(ID const& id) const {
     return (self->points_.find(id) != self->points_.end());
   }
-  
+
   /*!
    * @brief 回傳dimension
    */
   size_t dimension() const {
     return self->dimension_;
   }
-  
+
   /*!
    * @brief 設定dimension, 並且清空資料
    */
@@ -122,7 +118,7 @@ public:
     clear();
     return dimension();
   }
-  
+
   /*!
    * @brief 設定dimension, 並且針對每個identity point指定重設dimension
    */
@@ -134,14 +130,14 @@ public:
     }
     return dimension();
   }
-  
+
   /*!
    * @brief 取得所有identity points
    */
   IdentityPointsMap const& identityPoints() const {
     return self->points_;
   }
-  
+
   /*!
    * @brief 設定所有identity points
    */
@@ -149,7 +145,7 @@ public:
     clear();
     return identityPointsAdd(points);
   }
-  
+
   /*!
    * @brief 加入identity Points
    */
@@ -159,7 +155,7 @@ public:
     }
     return identityPoints();
   }
-  
+
   /*!
    * @brief 移除identity Points
    */
@@ -170,14 +166,14 @@ public:
     }
     return identityPoints();
   }
-  
+
   /*!
    * @brief 取得一個identity point
    */
   Vector<Scalar> identityPoint(ID const& id) const {
     return (exist(id) ? self->points_.find(id)->second : Vector<Scalar>());
   }
-  
+
   /*!
    * @brief 修改一個identity point
    */
@@ -187,7 +183,7 @@ public:
     }
     return identityPoint(id);
   }
-  
+
   /*!
    * @brief 新增一個identity point
    */
@@ -197,42 +193,99 @@ public:
     }
     return identityPoint(id);
   }
-  
+
   /*!
    * @brief 刪除一個identity point
    */
   void identityPointDel(ID const& id) {
     self()->points_.erase(id);
   }
-  
+
   /*!
    * @brief 取得一個identity point, non-constant reference
    */
   Vector<Scalar>& identityPointGet(ID const& id) {
     return self()->points_[id];
   }
-  
+
   /*!
    * @brief same as \c copyFrom(b)
    */
   IdentityPoints& operator=(IdentityPoints const& b) {
     return copyFrom(b);
   }
-  
+
   /*! @brief 將資料寫入檔案
    *
    *  @note 未完成
    */
   bool write(FILE* f, bool bin, unsigned int fg) const {
-    return false;
+    if (bin) {
+      long dim, ct;
+      if (fwrite(&(dim = dimension()), sizeof(dim), 1, f) < 1) return false;
+      if (fwrite(&(ct = size()), sizeof(ct), 1, f) < 1) return false;
+      for (IdentityPointsMapIterK
+           it = identityPoints().begin(), ed = identityPoints().end();
+           it != ed; ++it) {
+        double tmp;
+        if (fwrite(&(tmp = it->first), sizeof(tmp), 1, f) < 1) return false;
+        for (long i = 0; i < dim; ++i) {
+          if (fwrite(&(tmp = it->second(i)), sizeof(tmp), 1, f) < 1)
+            return false;
+        }
+      }
+    }
+    else {
+      if (fprintf(f, "%ld %lu\n", dimension(), size()) < 1) return false;
+      for (IdentityPointsMapIterK
+           it = identityPoints().begin(), ed = identityPoints().end();
+           it != ed; ++it) {
+        if (fprintf(f, "%f   ", (double)it->first) < 1) return false;
+        for (long i = 0, I = dimension(); i < I; ++i) {
+          if (fprintf(f, "%f ", (double)it->second(i)) < 1) return false;
+        }
+        fprintf(f, "\n");
+      }
+    }
+    return true;
   }
-  
+
   /*! @brief 將資料讀入
    *
    *  @note 未完成
    */
   bool read(FILE* f, bool bin, unsigned int fg) {
-    return false;
+    long dim, ct;
+    if (bin) {
+      if (fread(&dim, sizeof(dim), 1, f) < 1) return false;
+      dimension(dim);
+      if (fread(&ct, sizeof(ct), 1, f) < 1) return false;
+      double id, tt;
+      Vector<Scalar> tmp(dim, 0);
+      for (int i = 0; i < ct; ++i) {
+        if (fread(&id, sizeof(id), 1, f) < 1) return false;
+        for (size_t j = 0, J = dim; j < J; ++j) {
+          if (fread(&tt, sizeof(tt), 1, f) < 1) return false;
+          tmp.entry(j, tt);
+        }
+        identityPointAdd((ID)id, tmp);
+      }
+    }
+    else {
+      if (fscanf(f, "%ld %ld", &dim, &ct) < 2) return false;
+      dimension(dim);
+      double id, tt;
+      Vector<Scalar> tmp(dim, 0);
+      for (int i = 0; i < ct; ++i) {
+        if (fscanf(f, "%lf", &id) < 1) return false;
+        for (int j = 0, J = dim; j < J; ++j) {
+          if (fscanf(f, "%lf", &tt) < 1) return false;
+          tmp.entry(j, tt);
+        }
+        identityPointAdd((ID)id, tmp);
+      }
+    }
+    return true;
   }
 
   /*! @brief new一個自己
@@ -242,7 +295,7 @@ public:
   ObjBase* create() const {
     return new IdentityPoints();
   }
-  
+
   /*! @brief 複製資料
    *
    * 輸入型別是 \c ObjBase \c const*
@@ -255,16 +308,15 @@ public:
   ObjBase* copyFrom(ObjBase const* b) {
     return &(copyFrom(*(IdentityPoints*)b));
   }
-  
+
   /*! @brief 回傳class的type
    *
    *  @return \c char \c const\c * 形式的typename
    */
-  char const* ctype() const{
-    static char const* ptr = typeid(*this).name();
-    return ptr;
+  char const* ctype() const {
+    return typeid(*this).name();
   }
-  
+
   /*! @brief 回傳class的type
    *
    *  @return \c std::string 形式的typename
diff --git a/meowpp/gra/Photo.h b/meowpp/gra/Photo.h
index 15768aa..1837727 100644
--- a/meowpp/gra/Photo.h
+++ b/meowpp/gra/Photo.h
@@ -1,8 +1,11 @@
 #ifndef   gra_Photo_H__
 #define   gra_Photo_H__
 
+#include "Bitmap.h"
+
 #include "../Self.h"
 
+#include "../geo/Vectors.h"
 #include "../math/utility.h"
 #include "../math/Matrix.h"
 #include "../math/Transformations.h"
@@ -20,7 +23,7 @@ namespace meow {
 /*!
  * @brief 底片 
  * 
- * 基本上就是一個 \c Photo 加上 \c focal
+ * 基本上就是一個 \c Bitmap 加上 \c focal
  * 
  * @author cat_leopard
  */
@@ -33,13 +36,9 @@ private:
     PhotoProjection<double> proj_;
     Myself(): proj_(3) {
     }
-    ~Myself() {
+    Myself(Myself const& b): bmp_(b.bmp_), c_(b.c_), proj_(b.proj_) {
     }
-    Myself& copyFrom(Myself const& b) {
-      bmp_ .copyFrom(b. bmp_);
-      c_   .copyFrom(b.   c_);
-      proj_.copyFrom(b.proj_);
-      return *this;
+    ~Myself() {
     }
   };
   
@@ -57,7 +56,7 @@ public:
    * 
    * focal 預設為 1
    */
-  Photo(): self(true) {
+  Photo(): self() {
     self()->proj_.focal(1.0);
   }
   
@@ -68,8 +67,7 @@ public:
    * 
    * @param [in] b 資料來源
    */
-  Photo(Photo const& b): self(false) {
-    copyFrom(b);
+  Photo(Photo const& b): self(b.self, Self<Myself>::COPY_FROM) {
   }
   
   /*!
@@ -79,7 +77,7 @@ public:
    * 
    * @param [in] bmp 給定的圖片
    */
-  Photo(Bitmap<Pixel> const& bmp): self(true) {
+  Photo(Bitmap<Pixel> const& bmp): self() {
     reset(bmp);
   }
   
@@ -91,7 +89,7 @@ public:
    * @param [in] bmp 給定的圖片
    * @param [in] f   給定的焦距
    */
-  Photo(Bitmap<Pixel> const& bmp, double f): self(true) {
+  Photo(Bitmap<Pixel> const& bmp, double f): self() {
     reset(bmp, f);
   }
   
@@ -104,8 +102,7 @@ public:
    * @param [in] f   給定的焦距
    * @param [in] c   中心點作標
    */
-  Photo(Bitmap<Pixel> const& bmp, double f, Vector2D<double> const& c):
-  self(true) {
+  Photo(Bitmap<Pixel> const& bmp, double f, Vector2D<double> const& c): self() {
     reset(bmp, f, c);
   }
   
@@ -196,7 +193,7 @@ public:
    * @return 新的 \c bitmap
    */
   Bitmap<Pixel> const& bitmap(Bitmap<Pixel> const& bmp) {
-    self()->bmp_ = bmp;
+    self()->bmp_.copyFrom(bmp);
     return bitmap();
   }
   
@@ -219,6 +216,23 @@ public:
   }
   
   /*!
+   * @brief 回傳相應的 photo projection
+   */
+  PhotoProjection<double> projection() const {
+    return self->proj_;
+  }
+  
+  /*!
+   * @brief 設定 photo projection
+   */
+  PhotoProjection<double> projection(PhotoProjection<double> const& p) {
+    if (p.dimension() == 3) {
+      self()->proj_ = p;
+    }
+    return projection();
+  }
+  
+  /*!
    * @brief 取得照片中心點底片座標
    *
    * @return 一個二維vector
@@ -299,6 +313,7 @@ public:
    * @return \c true/false
    */
   bool inside(Vector3D<double> const& p) const {
+    if (p.z() > 0) return false;
     return inside(Vector2D<double>(self->proj_.transformate(p.matrix())));
   }
   
@@ -352,7 +367,18 @@ public:
    *  @note 未完成
    */
   bool write(FILE* f, bool bin, unsigned int fg) const {
-    return false;
+    if (bitmap().write(f, bin, fg) == false) return false;
+    if (bin) {
+      double tmp;
+      if (fwrite(&(tmp = center().x()), sizeof(tmp), 1, f) < 1) return false;
+      if (fwrite(&(tmp = center().y()), sizeof(tmp), 1, f) < 1) return false;
+      if (fwrite(&(tmp = focal()), sizeof(tmp), 1, f) < 1) return false;
+    }
+    else {
+      if (fprintf(f, "%f %f\n", center().x(), center().y()) < 2) return false;
+      if (fprintf(f, "%f\n", focal()) < 1) return false;
+    }
+    return true;
   }
   
   /*! @brief 將資料讀入
@@ -360,12 +386,23 @@ public:
    *  @note 未完成
    */
   bool read(FILE* f, bool bin, unsigned int fg) {
-    return false;
+    if (bitmapGet().read(f, bin, fg) == false) return false;
+    double tmp[3];
+    if (bin) {
+      if (fread(tmp, sizeof(double), 3, f) < 3) return false;
+    }
+    else {
+      if (fscanf(f, "%lf %lf %lf", tmp + 0, tmp + 1, tmp + 2) < 3) return false;
+    }
+    centerGet().x(tmp[0]);
+    centerGet().y(tmp[1]);
+    focal(tmp[2]);
+    return true;
   }
 
   /*! @brief new一個自己
    *
-   *  @return 一個new出來的Bitmap<Pixel>
+   *  @return 一個new出來的Photo<Pixel>
    */
   ObjBase* create() const {
     return new Photo();
@@ -389,8 +426,7 @@ public:
    *  @return \c char \c const\c * 形式的typename
    */
   char const* ctype() const{
-    static char const* ptr = typeid(*this).name();
-    return ptr;
+    return typeid(*this).name();
   }
   
   /*! @brief 回傳class的type
@@ -402,6 +438,6 @@ public:
   }
 };
   
-}
+} // meow
 
 #endif // gra_Photo_H__
diff --git a/meowpp/gra/WatchBall.h b/meowpp/gra/WatchBall.h
index c72e2a5..f3d2044 100644
--- a/meowpp/gra/WatchBall.h
+++ b/meowpp/gra/WatchBall.h
@@ -43,14 +43,13 @@ public:
   /*!
    * @brief constructor
    */
-  WatchBall(): self(true) {
+  WatchBall(): self() {
   }
   
   /*!
    * @brief copy constructor
    */
-  WatchBall(WatchBall const& b): self(false) {
-    copyFrom(b);
+  WatchBall(WatchBall const& b): self(b.self, Self<Myself>::COPY_FROM) {
   }
   
   /*!
@@ -176,12 +175,12 @@ public:
     Bitmap<Pixel> ret(height, width, Pixel(0));
     for (size_t i = 0; i < height; ++i) {
       for (size_t j = 0; j < width; ++j) {
-        double theta = (1.0 * j /  width - 0.5) * 2 * PI;
+        double theta = (1.0 * j / width - 0.5) * 2 * PI;
         double phi   = asin(-(1.0 * i / height - 0.5) * 2.0);
         ret.pixel(i, j, color(Vector3D<double>(
-              sin(-theta) * cos(phi),
+              sin(theta) * cos(phi),
               sin(phi),
-              cos(-theta) * cos(phi)
+              -cos(theta) * cos(phi)
         )));
       }
     }
diff --git a/meowpp/math/LinearTransformation.h b/meowpp/math/LinearTransformation.h
index 86724b7..f51630e 100644
--- a/meowpp/math/LinearTransformation.h
+++ b/meowpp/math/LinearTransformation.h
@@ -71,15 +71,13 @@ protected:
   }
   
   /*!
-   * @brief Setup the matrix
-   *
-   * @param [in] m matrix
-   * @return new matrix
+   * @brief setup the matrix
    */
-  Matrix<Scalar> const& matrix(Matrix<Scalar> const& m) {
+  virtual Matrix<Scalar> const& matrix(Matrix<Scalar> const& m) {
     matrix_.copyFrom(m);
-    return matrix_;
+    return matrix();
   }
+  
 public:
   /*!
    * Destructor
@@ -106,6 +104,7 @@ public:
   }
 };
 
-}
+
+} // meow
 
 #endif // math_LinearTransformation_H__
diff --git a/meowpp/math/LinearTransformations.h b/meowpp/math/LinearTransformations.h
index c7295ab..e882a6c 100644
--- a/meowpp/math/LinearTransformations.h
+++ b/meowpp/math/LinearTransformations.h
@@ -23,37 +23,41 @@ private:
     Vector3D<Scalar> theta_;
     bool              need_;
     
-    Myself() {
+    Myself(): theta_(0, 0, 0), need_(true) {
     }
-    ~Myself() {
+    
+    Myself(Myself const& b): theta_(b.theta_), need_(b.need_) {
     }
-    Myself& copyFrom(Myself const& b) {
-      theta_ = b.theta_;
-      need_  = b.need_;
-      return *this;
+    
+    ~Myself() {
     }
   };
   
-  Self<Myself> const& self;
+  Self<Myself> const self;
   
   void calcMatrix() const {
     if (self->need_) {
-      Vector3D<double> axis (self->theta_.normalize());
-      double           angle(self->theta_.length());
-      double           cs(cos(angle / 2.0));
-      double           sn(sin(angle / 2.0));
+      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));
 
-      Matrix<Scalar> tmp(3, 3, Scalar(0.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)->matrix(tmp);
+        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;
     }
   }
@@ -62,21 +66,14 @@ public:
   /*!
    * Constructor with no rotation
    */
-  Rotation3D(): LinearTransformation<Scalar>(3u, 3u, 3u),
-  self(true) {
-    self()->theta_.x(Scalar(0));
-    self()->theta_.y(Scalar(0));
-    self()->theta_.z(Scalar(0));
-    self()->need_ = true;
-    calcMatrix();
+  Rotation3D(): LinearTransformation<Scalar>(3u, 3u, 3u), self() {
   }
   
   /*!
    * Constructor and copy data
    */
   Rotation3D(Rotation3D const& b): LinearTransformation<Scalar>(b),
-  self(false) {
-    copyFrom(b);
+  self(b.self, Self<Myself>::COPY_FROM) {
   }
   
   /*!
@@ -214,7 +211,7 @@ public:
    * @return the output matrix
    */
   Matrix<Scalar> transformate(Matrix<Scalar> const& x) const {
-    if (self->need_) calcMatrix();
+    calcMatrix();
     return LinearTransformation<Scalar>::matrix() * x;
   }
   
@@ -244,7 +241,7 @@ public:
    * @return a matrix
    */
   Matrix<Scalar> jacobian(Matrix<Scalar> const& x) const {
-    if (self->need_) calcMatrix();
+    calcMatrix();
     return LinearTransformation<Scalar>::matrix();
   }
   
@@ -321,7 +318,7 @@ public:
    * @return a matrix
    */
   Matrix<Scalar> jacobian(Matrix<Scalar> const& x, size_t i) const {
-    if (self->need_) calcMatrix();
+    calcMatrix();
     Matrix<Scalar> mid(3u, 3u, Scalar(0.0));
     if (i == 0) {
       mid.entry(1, 2, Scalar(-1.0));
@@ -345,8 +342,7 @@ public:
    * @return the output vector
    */
   Matrix<Scalar> transformateInv(Matrix<Scalar> const& x) const {
-    if (self->need_) calcMatrix();
-    return LinearTransformation<Scalar>::matrix().transpose() * x;
+    return matrixInv() * x;
   }
   
   /*!
@@ -356,8 +352,7 @@ public:
    * @return a matrix
    */
   Matrix<Scalar> jacobianInv(Matrix<Scalar> const& x) const {
-    if (self->need_) calcMatrix();
-    return LinearTransformation<Scalar>::matrix().transpose();
+    return matrixInv();
   }
   
   /*!
@@ -368,7 +363,7 @@ public:
    * @return a matrix
    */
   Matrix<Scalar> jacobianInv(Matrix<Scalar> const& x, size_t i) const {
-    if (self->need_) calcMatrix();
+    calcMatrix();
     Matrix<Scalar> mid(3u, 3u, Scalar(0.0));
     if (i == 0) {
       mid.entry(1, 2, Scalar(-1.0));
@@ -382,7 +377,8 @@ public:
       mid.entry(0, 1, Scalar(-1.0));
       mid.entry(1, 0, Scalar( 1.0));
     }
-    return mid.transpose() * matrixInv() * x;
+    return matrixInv() * mid.transpose() * x;
+    return (-mid) * matrixInv() * x;
   }
 
   /*!
@@ -393,7 +389,7 @@ public:
    * @return a matrix
    */
   Matrix<Scalar> matrixInv() const {
-    if (self->need_) calcMatrix();
+    calcMatrix();
     return LinearTransformation<Scalar>::matrix().transpose();
   }
   
diff --git a/meowpp/math/Matrix.h b/meowpp/math/Matrix.h
index b4b4853..09832d8 100644
--- a/meowpp/math/Matrix.h
+++ b/meowpp/math/Matrix.h
@@ -1,8 +1,6 @@
 #ifndef   math_Matrix_H__
 #define   math_Matrix_H__
 
-#include "utility.h"
-
 #include "../Self.h"
 
 #include <vector>
@@ -18,25 +16,32 @@ namespace meow {
  */
 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():
+    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;
     }
-    Myself& copyFrom(Myself const& m) {
-      rows_    = m.   rows_;
-      cols_    = m.   cols_;
-      entries_ = m.entries_;
-      return *this;
-    }
   };
+
   Self<Myself> const self;
 public:
   /*!
@@ -45,7 +50,7 @@ public:
    * Create an empty matrix with size \b 0x0.
    * In other world, create an \b invalid matrix
    */
-  Matrix(): self(true) { }
+  Matrix(): self() { }
 
   /*!
    * @brief constructor
@@ -54,7 +59,8 @@ public:
    *
    * @param [in] m another matrix
    */
-  Matrix(Matrix const& m): self(false) { self().copyFrom(m.self); }
+  Matrix(Matrix const& m): self(m.self, Self<Myself>::COPY_FROM) {
+  }
 
   /*!
    * @brief constructor
@@ -65,7 +71,8 @@ public:
    * @param [in] c number of columns
    * @param [in] e inital entry
    */
-  Matrix(size_t r, size_t c, Entry const& e): self(true) { reset(r, c, e); }
+  Matrix(size_t r, size_t c, Entry const& e): self(Myself(r, c, e)) {
+  }
 
   //! @brief destructor
   ~Matrix() { }
@@ -152,8 +159,7 @@ public:
    */
   size_t cols(size_t c, Entry const& e) {
     if (c != cols()) {
-      Self<Myself> const old(false);
-      old().copyFrom(self);
+      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++) {
@@ -169,7 +175,7 @@ public:
 
   /*!
    * @brief resize
-   * 
+   *
    * Resize to \a r x \a c, with new created entry be \a e
    *
    * @param [in] r number of rows
@@ -193,14 +199,19 @@ public:
     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] rFirst
+   * @param [in] rLast
    * @param [in] cFirst
-   * @param [in] cLast 
+   * @param [in] cLast
    * @param [in] e value
    * @return void
    */
@@ -218,7 +229,7 @@ public:
    * @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
@@ -228,7 +239,7 @@ public:
   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) {
+    if (rFirst == 0 && cFirst == 0) {
       Matrix ret(*this);
       ret.size(rLast + 1, cLast + 1, Entry(0));
       return ret;
@@ -340,6 +351,28 @@ public:
         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)
@@ -372,7 +405,7 @@ public:
     tmp.size(cols(), rows(), Entry(0));
     return tmp;
   }
-  
+
   //! @brief let itself become itself's inverse matrix
   Matrix& inversed() {
     copyFrom(inverse());
@@ -380,14 +413,14 @@ public:
   }
 
   //! @brief return itself's transpose matrix
-  Matrix  transpose () const {
+  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());
@@ -400,7 +433,7 @@ public:
     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++) {
@@ -432,53 +465,53 @@ public:
   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
index 0d9ee05..011db8c 100644
--- a/meowpp/math/Transformation.h
+++ b/meowpp/math/Transformation.h
@@ -4,6 +4,7 @@
 #include "Matrix.h"
 #include "../Self.h"
 
+#include <list>
 #include <cstdlib>
 
 namespace meow {
@@ -13,7 +14,7 @@ namespace meow {
  *
  * We define that the input and output form of our transformations all be
  * \b matrix . Some advance methods such as calculating jacobian matrix
- * will order that the input form must be a vector.
+ * will require that the input form must be a vector.
  * @author cat_leopard
  */
 template<class Scalar>
@@ -25,17 +26,22 @@ private:
     size_t outputRows_;
     size_t outputCols_;
     size_t      psize_;
-    
-    Myself& copyFrom(Myself const& b) {
-      inputRows_  = b. inputRows_;
-      inputCols_  = b. inputCols_;
-      outputRows_ = b.outputRows_;
-      outputCols_ = b.outputCols_;
-      psize_      = b.psize_;
-      return *this;
+
+    Myself() {
+    }
+    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:
   /*!
@@ -48,22 +54,18 @@ protected:
    */
   Transformation(size_t  inputRows, size_t  inputCols,
                  size_t outputRows, size_t outputCols,
-                 size_t psize): self(true) {
-    self()-> inputRows_ =  inputRows;
-    self()-> inputCols_ =  inputCols;
-    self()->outputRows_ = outputRows;
-    self()->outputCols_ = outputCols;
-    self()->psize_      = psize;
+                 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(false) {
-    copyFrom(b);
+  Transformation(Transformation const& b):
+  self(b.self, Self<Myself>::COPY_FROM) {
   }
-  
+
   /*!
    * @brief Copy from the specified one
    *
@@ -74,9 +76,9 @@ protected:
     self().copyFrom(b.self);
     return *this;
   }
-  
+
   /*!
-   * @brief Ceference from the specified one
+   * @brief reference from the specified one
    *
    * @param [in] b The specified one
    * @return \c *this
@@ -127,7 +129,7 @@ public:
   size_t outputCols() const {
     return self->outputCols_;
   }
-  
+
   /*!
    * @brief Return the number of parameters.
    *
@@ -230,6 +232,6 @@ public:
   }
 };
 
-}
+} // meow
 
 #endif // math_Transformation_H__
diff --git a/meowpp/math/Transformations.h b/meowpp/math/Transformations.h
index b712558..b9583f9 100644
--- a/meowpp/math/Transformations.h
+++ b/meowpp/math/Transformations.h
@@ -50,17 +50,14 @@ template<class Scalar>
 class BallProjection: public Transformation<Scalar> {
 private:
   struct Myself {
-    Scalar    radius_;
     size_t dimension_;
+    Scalar    radius_;
 
-    Myself() {
+    Myself(size_t d): dimension_(1), radius_(1) {
     }
-    ~Myself() {
+    Myself(size_t d, Scalar const& r): dimension_(d), radius_(r) {
     }
-    Myself& copyFrom(Myself const& b) {
-      radius_    = b.radius_;
-      dimension_ = b.dimension_;
-      return *this;
+    Myself(Myself const& m): dimension_(m.dimension_), radius_(m.radius_) {
     }
   };
   
@@ -70,18 +67,16 @@ public:
    * Constructor, copy settings from given BallProjection
    * @param [in] b another ball projection class
    */
-  BallProjection(BallProjection const& b):
-  Transformation<Scalar>(b), self(false) {
-    copyFrom(b);
+  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): self(true),
-  Transformation<Scalar>(d, 1, d, 1, 1) {
-    self()->dimension_ = d;
+  BallProjection(size_t d): Transformation<Scalar>(d, 1, d, 1, 1),
+  self(Myself(d)) {
     radius(1);
   }
   
@@ -90,9 +85,8 @@ public:
    * @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(true) {
-    self()->dimension_ = d;
+  BallProjection(size_t d, Scalar const& r): Transformation<Scalar>(d,1,d,1,1),
+  self(Myself(d, r)) {
     radius(r);
   }
   
@@ -334,41 +328,40 @@ private:
     
     Myself() {
     }
-    ~Myself() {
+    
+    Myself(size_t d, Scalar f): focal_(f), dimension_(d) {
     }
-    Myself& copyFrom(Myself const& b) {
-      focal_     = b.focal_;
-      dimension_ = b.dimension_;
-      return *this;
+    
+    Myself(Myself const& b): focal_(b.focal_), dimension_(b.dimension_) {
+    }
+    
+    ~Myself() {
     }
   };
   
-  Self<Myself> const& self;
+  Self<Myself> const self;
 public:
   /*!
    * Constructor, focal = 1
    */
   PhotoProjection(size_t dimension): 
-  Transformation<Scalar>(dimension, 1, dimension, 1, 1), self(true) {
-    self()->dimension_ = dimension;
-    focal(1);
+  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(true) {
-    self()->dimension_ = dimension;
-    focal(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(false) {
-    self().copyFrom(p.self);
+  PhotoProjection(PhotoProjection const& p): Transformation<Scalar>(p),
+  self(p.self, Self<Myself>::COPY_FROM) {
   }
   
   /*!
@@ -510,7 +503,6 @@ public:
    *  .
    * Then the jacobian matrix is like below: \n
    * \f[
-   *   f \times
    *   \left[
    *     \begin{array}{c}
    *       \frac{-x_1}{x_N} \\
diff --git a/meowpp/math/Vector.h b/meowpp/math/Vector.h
index d387c2b..caa64fd 100644
--- a/meowpp/math/Vector.h
+++ b/meowpp/math/Vector.h
@@ -3,7 +3,6 @@
 
 #include "../Self.h"
 #include "Matrix.h"
-#include "utility.h"
 
 #include <vector>
 
@@ -18,6 +17,9 @@ namespace meow {
  */
 template<class Scalar>
 class Vector {
+public:
+  typedef typename Matrix<Scalar>::EntryRefK ScalarRefK;
+  typedef typename Matrix<Scalar>::EntryRef  ScalarRef ;
 private:
   Matrix<Scalar> matrix_;
 public:
@@ -26,9 +28,9 @@ public:
    *
    * With \b dimension=0, which means \b invalid.
    */
-  Vector(){
+  Vector() {
   }
-  
+
   /*!
    * @brief constructor
    *
@@ -36,10 +38,9 @@ public:
    *
    * @param [in] v another vector
    */
-  Vector(Vector const& v) {
-    matrix_.copyFrom(v.matrix_);
+  Vector(Vector const& v): matrix_(v.matrix_) {
   }
-  
+
   /*!
    * @brief constructor
    *
@@ -47,22 +48,9 @@ public:
    *
    * @param [in] m matrix
    */
-  Vector(Matrix<Scalar> const& m) {
-    matrix_.copyFrom(m.col(0));
-  }
-  
-  /*!
-   * @brief constructor
-   *
-   * From matrix's \a i-th column
-   *
-   * @param [in] m matrix
-   * @param [in] i i-th
-   */
-  Vector(Matrix<Scalar> const& m, size_t i) {
-    matrix_.copyFrom(m.col(i));
+  Vector(Matrix<Scalar> const& m): matrix_(m.col(0)) {
   }
-  
+
   /*!
    * @brief constructor
    *
@@ -70,27 +58,25 @@ public:
    *
    * @param [in] v vector
    */
-  Vector(std::vector<Scalar> const& v) {
-    matrix_.size(v.size(), 1, Scalar(0));
+  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_.reset(d, 1, e);
+  Vector(size_t d, Scalar const& e): matrix_(d, 1, e) {
   }
 
   //! @brief destructor
-  ~Vector(){
+  ~Vector() {
   }
 
   //! @brief copy from ...
@@ -114,7 +100,7 @@ public:
   size_t dimension() const {
     return matrix_.rows();
   }
-  
+
   /*!
    * @brief resize the dimension
    *
@@ -126,7 +112,7 @@ public:
     matrix_.rows(d, s);
     return dimension();
   }
-  
+
   /*!
    * @brief Return whether \c dimension>0 is true or not
    * @return \c true/false
@@ -134,12 +120,12 @@ public:
   bool valid() const {
     return (dimension() > 0);
   }
-  
+
   //! @brief return \a i -th entry
   Scalar entry(size_t i) const {
     return matrix_.entry(i, 0);
   }
-  
+
   /*!
    * @brief change \a i -th entry
    *
@@ -151,6 +137,11 @@ public:
     return entry(i);
   }
   
+  //! @brief return \a i -th entry with non-constant type
+  ScalarRef entryGet(size_t i) {
+    return matrix_.entryGet(i);
+  }
+
   /*!
    * @brief change \a i -th to \a j -th entries
    *
@@ -159,11 +150,11 @@ public:
    * @param [in] s new value
    */
   void entries(size_t i, size_t j, Scalar const& s) {
-    for (size_t it = i; it <= j; it++) {
+    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));
@@ -173,7 +164,7 @@ public:
   Vector positive() const {
     return *this;
   }
-  
+
   //! @brief return -\a (*this)
   Vector negative() const {
     return Vector(matrix_.negative());
@@ -183,7 +174,7 @@ public:
   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_));
@@ -193,7 +184,7 @@ public:
   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));
@@ -208,17 +199,17 @@ public:
   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());
@@ -229,48 +220,48 @@ public:
   Vector& operator=(Vector const& v) {
     return copyFrom(v);
   }
-  
+
   //! @brief same as entry(i)
   Scalar operator()(size_t i) const {
     return entry(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
index 4a47d33..5b67594 100644
--- a/meowpp/math/methods.h
+++ b/meowpp/math/methods.h
@@ -92,7 +92,7 @@ inline std::vector<Data> ransac(std::vector<Data> const& data,
 }
 
 
-/*!
+/*
  * @brief Run the \b Levenberg-Marquardt method to solve a non-linear 
  *        least squares problem.
  *
@@ -159,23 +159,72 @@ inline std::vector<Data> ransac(std::vector<Data> const& data,
  *
  * @author cat_leopard
  */
-template<class Scalar, class F, class J, class I, class Stop>
-inline Vector<Scalar> levenbergMarquardt(F const& func,
-                                         J const& jaco,
-                                         I const& iden,
+template<class Scalar, class Function>
+inline Vector<Scalar> levenbergMarquardt(Function       const& f,
                                          Vector<Scalar> const& init,
-                                         Stop           const& stop,
-                                         int counter = -1) {
-  Vector<Scalar> ans(init);
-  for(Vector<Scalar> rv;
-      !stop((rv = func(ans)).length2()) && counter != 0; counter--) {
-    Matrix<Scalar> j(jaco(ans)), jt(j.transpose());
-    Matrix<Scalar> i(iden(ans));
-    ans = ans - Vector<Scalar>((jt * j + i).inverse() * jt * rv.matrix());
+                                         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
index 43b4785..d153b68 100644
--- a/meowpp/math/utility.h
+++ b/meowpp/math/utility.h
@@ -6,7 +6,7 @@
 #include <algorithm>
 #include <cmath>
 
-namespace meow{
+namespace meow {
 
 //! 圓周率...
 static const double PI = 3.14159265358979323846264338327950288;
@@ -15,7 +15,7 @@ static const double PI = 3.14159265358979323846264338327950288;
  * @brief 如果abs(輸入的數值) < eps, 則回傳0, 否則回傳輸入的數值
  */
 template<class T>
-inline T noEPS(T value, T eps = 1e-9){
+inline T noEPS(T value, T eps = 1e-9) {
   T epsp((eps < T(0)) ? -eps : eps);
   return ((value < -epsp || value > epsp) ? value : T(0));
 }
@@ -24,7 +24,7 @@ inline T noEPS(T value, T eps = 1e-9){
  * @brief \c (value-lower)/(upper-lower)
  */
 template<class T>
-inline T normalize(T lower, T upper, T value){
+inline T normalize(T lower, T upper, T value) {
   return (value - lower) / (upper - lower);
 }
 
@@ -32,7 +32,7 @@ inline T normalize(T lower, T upper, T value){
  * @brief \c (lower+_ratio*(upper-lower))
  */
 template<class T>
-inline T denormalize(T lower, T upper, T _ratio){
+inline T denormalize(T lower, T upper, T _ratio) {
   return lower + _ratio * (upper - lower);
 }
 
@@ -40,7 +40,7 @@ inline T denormalize(T lower, T upper, T _ratio){
  * @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){
+inline T ratioMapping(T l1, T u1, T m1, T l2, T u2) {
   return denormalize(l2, u2, normalize(l1, u1, m1));
 }
 
@@ -48,15 +48,23 @@ inline T ratioMapping(T l1, T u1, T m1, T l2, T u2){
  * @brief \c std::min(mx,std::max(mn,v))
  */
 template<class T>
-inline T inRange(T const& mn, T const& mx, T const& v){
+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){
+inline T squ(T const& x) {
   return x * x;
 }
 
@@ -64,7 +72,7 @@ inline T squ(T const& x){
  * @brief \c x*x*x
  */
 template<class T>
-inline T cub(T const& x){
+inline T cub(T const& x) {
   return x * x * x;
 }
 
@@ -72,22 +80,22 @@ inline T cub(T const& x){
  * @brief 只將 \c sigs 個標準差以內的數據拿來取平均
  */
 template<class T>
-inline double average(T const& beg, T const& end, double sigs){
+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++){
+  for (T it = beg; it != end; ++it, ++N) {
     av += *it;
   }
   av /= N;
   double sig = 0;
-  for(T it = beg; it != end; it++){
+  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){
+  for (T it = beg; it != end; ++it) {
+    if (lower <= *it && *it <= upper) {
       ret += *it;
       retn++;
     }
@@ -99,30 +107,30 @@ inline double average(T const& beg, T const& end, double sigs){
  * @brief 只將 \c sigs 個標準差以內的數據拿來取平均, 不過這次用 \c p 來加權平均
  */
 template<class T>
-inline double average(T const& beg, T const& end, T const& p, double sigs){
+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++){
+  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++){
+  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++){
+  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){
+  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;
+  if (retn <= 1e-10) return av;
   return ret / retn;
 }
 
@@ -130,10 +138,10 @@ inline double average(T const& beg, T const& end, T const& p, double sigs){
  * @brief 就只是個取絕對值
  */
 template<class T>
-inline T tAbs(T const& t){
+inline T tAbs(T const& t) {
   return (t < 0 ? -t : t);
 }
 
-}
+} // meow
 
 #endif // math_utility_H__
diff --git a/meowpp/oo/ObjArray.h b/meowpp/oo/ObjArray.h
index 4417524..b58b89c 100644
--- a/meowpp/oo/ObjArray.h
+++ b/meowpp/oo/ObjArray.h
@@ -37,19 +37,19 @@ private:
 public:
   ObjArray(): self(true) {
   }
-  
+
   ObjArray(ObjArray const& a): self(false) {
     self().copyFrom(a.self);
   }
-  
+
   ObjArray(std::vector<T> const& a): self(true) {
     self()->array_ = a;
   }
-  
+
   ObjArray(size_t sz, T const& e): self(true) {
     self()->array_.resize(sz, e);
   }
-  
+
   ~ObjArray() {
   }
 
@@ -79,7 +79,7 @@ public:
     self()->array_.resize(res);
     return size();
   }
-  
+
   void clear() {
     self()->array_.clear();
   }
@@ -96,7 +96,7 @@ public:
     self()->array_.push_back(e);
     return entry(size() - 1);
   }
-  
+
   bool popBack() {
     if (empty()) return false;
     self()->array_.pop_back();
@@ -106,11 +106,11 @@ public:
   ObjArray& operator=(ObjArray const& a) {
     return copyFrom(a);
   }
-  
+
   T const& operator[](size_t i) const {
     return self->array_[i];
   }
-  
+
   T& operator[](size_t i) {
     return self()->array_[i];
   }
@@ -128,7 +128,7 @@ public:
     }
     return true;
   }
-  
+
   bool read(FILE* f, bool bin, unsigned int fg) {
     size_t sz;
     if (bin) {
@@ -143,19 +143,19 @@ public:
     }
     return true;
   }
-  
+
   ObjBase* create() const {
     return new ObjArray();
   }
-  
+
   ObjBase* copyFrom(ObjBase const* b) {
     return &(copyFrom(*b));
   }
-  
+
   char const* ctype() const {
     return typeid(*this).name();
   }
-  
+
   std::string type() const {
     return std::string(ctype());
   }
diff --git a/meowpp/oo/ObjBase.h b/meowpp/oo/ObjBase.h
index 974fdf1..0ae9427 100644
--- a/meowpp/oo/ObjBase.h
+++ b/meowpp/oo/ObjBase.h
@@ -2,11 +2,10 @@
 #define   oo_ObjBase_H__
 
 #include <cstdio>
-
 #include <typeinfo>
 #include <string>
 
-namespace meow{
+namespace meow {
 
 /*!
  * @brief 一切物件的Base, 並要求每個物件都要有read, write, create, ... 等功能
@@ -18,26 +17,80 @@ protected:
   ObjBase(){ }
 public:
   virtual ~ObjBase(){ }
-  //
-  virtual bool write(FILE* f,bool bin,unsigned int fg) const { return false; }
-  virtual bool  read(FILE* f,bool bin,unsigned int fg)       { return false; }
-  //
-  virtual ObjBase* create() const { return NULL; }
-  virtual ObjBase* copyFrom(ObjBase const* b) { (*this) = (*b); return this; }
-  //
-  virtual char const* ctype() const{
-    static char const* ptr = typeid(*this).name();
-    return ptr;
-  }
-  virtual std::string type() const{ return std::string(ctype()); }
-  //
-  static char const* ctypeBase(){
-    static char const* ptr = typeid(ObjBase).name();
-    return ptr;
-  }
-  static std::string typeBase(){ return std::string(ctypeBase()); }
+  
+  /*!
+   * @brief 將物件寫入檔案, 預設implement為直接回傳 \c false
+   *
+   * @param [in] f   檔案
+   * @param [in] bin 是否為binary模式
+   * @param [in] fg  使用者自訂的argument
+   * @return 成功或失敗
+   */
+  virtual bool write(FILE* f, bool bin, unsigned int fg) const {
+    return false;
+  }
+  
+  /*!
+   * @brief 將物件從檔案讀出, 預設implement為直接回傳 \c false
+   *
+   * @param [in] f   檔案
+   * @param [in] bin 是否為binary模式
+   * @param [in] fg  使用者自訂的argument
+   * @return 成功或失敗
+   */
+  virtual bool read(FILE* f, bool bin, unsigned int fg) {
+    return false;
+  }
+  
+  /*!
+   * @brief 回傳一個new出來的物件, 預設implement為直接回傳 \c NULL
+   */
+  virtual ObjBase* create() const {
+    return NULL;
+  }
+  
+  /*!
+   * @brief 複製, 預設使用operator=
+   *
+   * @param [in] b 資料來源
+   * @return \c this
+   */
+  virtual ObjBase* copyFrom(ObjBase const* b) {
+    (*this) = (*b);
+    return this;
+  }
+  
+  /*!
+   * @brief 用C-style string回傳這個class的type name
+   */
+  virtual char const* ctype() const {
+    return typeid(*this).name();
+  }
+  
+  /*!
+   * @brief 用std::string回傳這個class的type name
+   */
+  virtual std::string type() const {
+    static std::string s(ctype());
+    return s;
+  }
+  
+  /*!
+   * @brief 用C-style string回傳base的type name
+   */
+  static char const* ctypeBase() {
+    return typeid(ObjBase).name();
+  }
+  
+  /*!
+   * @brief 用std::string回傳base的type name
+   */
+  static std::string typeBase() {
+    static std::string s(ctypeBase());
+    return s;
+  }
 };
 
-}
+} // meow
 
 #endif // oo_ObjBase_H__
diff --git a/meowpp/oo/ObjSelector.h b/meowpp/oo/ObjSelector.h
index 7664a11..58d86a0 100644
--- a/meowpp/oo/ObjSelector.h
+++ b/meowpp/oo/ObjSelector.h
@@ -11,7 +11,7 @@
 #include <cstdlib>
 #include <cstdio>
 
-namespace meow{
+namespace meow {
 
 /*!
  * @brief 利用register的概念, 達到runtime用string選擇要new的class
@@ -186,7 +186,7 @@ public:
    * @brief 從檔案中讀取一個物件(該物件必須要有註冊過)
    */
   static ObjBase* read(FILE* f, bool binary) {
-    char name[2048];
+    static char name[2048];
     size_t len;
     unsigned int fg;
     if (binary) {
diff --git a/meowpp/oo/ObjTypes.h b/meowpp/oo/ObjTypes.h
index 4092d63..48a4a07 100644
--- a/meowpp/oo/ObjTypes.h
+++ b/meowpp/oo/ObjTypes.h
@@ -21,6 +21,8 @@ private:
     Type data_;
     Myself() {
     }
+    Myself(Type const& t): data_(t) {
+    }
     ~Myself() {
     }
     Myself copyFrom(Myself const& b) {
@@ -30,17 +32,15 @@ private:
   Self<data_> const self;
 public:
   //! @brief constructor
-  ObjType(): self(true) {
+  ObjType(): self() {
   }
   
-  //! @brief constructor, 並且copy資料
-  ObjType(ObjType const& a): self(false) {
-    self().copyFrom(a.self);
+  //! @brief constructor, 並且給值
+  ObyType(Type const& t): self(Myself(t)) {
   }
   
-  //! @brief constructor, 並且給值
-  ObyType(Type const& t): self(true) {
-    self().data_ = t;
+  //! @brief constructor, 並且copy資料
+  ObjType(ObjType const& a): self(a.self, COPY_FROM) {
   }
   
   ~ObjType() {
@@ -93,8 +93,7 @@ public:
   }
   
   char const* ctype() const {
-    static char const* ptr = typeid(*this).name();
-    return ptr;
+    return typeid(*this).name();
   }
   
   std::string type() const {
@@ -197,6 +196,6 @@ typedef ObjType<size_t     , ReaderWriter_size_t> ObjSizeT;
 typedef ObjType<double     , ReaderWriter_double> ObjDouble;
 typedef ObjType<std::string, ReaderWriter_string> ObjString;
 
-}
+} // meow
 
 #endif // oo_ObjType_H__
diff --git a/meowpp/utility.h b/meowpp/utility.h
index b430dd8..83c7554 100644
--- a/meowpp/utility.h
+++ b/meowpp/utility.h
@@ -11,13 +11,35 @@
 namespace meow {
 
 /*!
+ * @brief 有.from.first, .from.second, .to.first, .to.second
+ *
+ * @author cathook
+ */
+template<class F1, class F2 = F1, class T1 = F1, class T2 = T1>
+struct PairToPair {
+  std::pair<F1, F2> from;
+  std::pair<T1, T2> to;
+  
+  PairToPair() {
+  }
+  PairToPair(PairToPair const& pp): from(pp.from), to(pp.to) {
+  }
+  PairToPair(F1 const& f1, F2 const& f2, T1 const& t1, T2 const& t2):
+  from(f1, f2), to(t1, t2) {
+  }
+  bool operator==(PairToPair const& p) const {
+    return (from == p.from && to == p.to);
+  }
+};
+
+/*!
  * @brief 類似C的printf, 不過是將格式化的字串丟到 \c std::string 裡回傳
- * 
+ *
  * @param [in] fmt,... 同printf
  * @return 一個 \c std::string
  * @warning 目前格式化字串最長只支援8191個字元
  */
-inline std::string stringPrintf(char const * fmt, ...){
+inline std::string stringPrintf(char const * fmt, ...) {
   char str[8192];
   va_list args;
   va_start(args, fmt);
@@ -28,7 +50,7 @@ inline std::string stringPrintf(char const * fmt, ...){
 
 /*!
  * @brief 將輸入字串中的某個pattern取代成另一個pattern
- * 
+ *
  * @param [in] str  輸入的字串
  * @param [in] from 要被取代的pattern
  * @param [in] to   將要取代的pattern
@@ -37,10 +59,10 @@ inline std::string stringPrintf(char const * fmt, ...){
  */
 inline std::string stringReplace(std::string str,
                                  std::string const& from,
-                                 std::string const& to){
+                                 std::string const& to) {
   std::string out = str;
   int len = from.length();
-  for(size_t pos; (pos = out.find(from)) != std::string::npos; ){
+  for (size_t pos; (pos = out.find(from)) != std::string::npos; ) {
     out.replace(pos, len, to);
   }
   return out;
@@ -56,14 +78,14 @@ inline std::string stringReplace(std::string str,
  * @note 參數中所有的字串都是用\c cstring 來傳遞,
  *       也就是 \c char \c const* 型態
  */
-inline bool cstringEndWith(char const* str, int n, ...){
+inline bool cstringEndWith(char const* str, int n, ...) {
   int len = strlen(str);
   va_list args;
   va_start(args, n);
-  for(int i = 0; i < n; i++){
+  for (int i = 0; i < n; i++) {
     char const* arg = va_arg(args, char const*);
     int arglen = strlen(arg);
-    if(arglen <= len && strcmp(str + len - arglen, arg) == 0){
+    if (arglen <= len && strcmp(str + len - arglen, arg) == 0) {
       return true;
     }
   }
@@ -87,7 +109,7 @@ debugPrintf_(\
 inline void debugPrintf_(char const* file,
                          char const* func,
                          size_t      line,
-                         char const* msg){
+                         char const* msg) {
 #ifdef    DEBUG
   fprintf(stderr, "%s[%d] %s >> %s", file, line, func, msg);
 #endif // DEBUG
@@ -96,7 +118,7 @@ inline void debugPrintf_(char const* file,
 /*!
  * @brief 階層式輸出
  *
- * 像是printf, 但多了一個 \b 巢狀 的概念, 例如: 
+ * 像是printf, 但多了一個 \b 巢狀 的概念, 例如:
  *
  * @code
  * message1(level = 0)
@@ -120,7 +142,7 @@ inline void debugPrintf_(char const* file,
  * @param [in] fmt, ... 跟printf一樣
  * @return 無
  */
-inline void messagePrintf(int level_change, char const* fmt, ...){
+inline void messagePrintf(int level_change, char const* fmt, ...) {
   static int      level =  0;
   static int last_level = -5;
   char str[8192];
@@ -128,14 +150,15 @@ inline void messagePrintf(int level_change, char const* fmt, ...){
   va_start(args, fmt);
   vsnprintf(str, 8192, fmt, args);
   va_end(args);
-  if(last_level == 1 && level_change == -1){
+  if (last_level == 1 && level_change == -1) {
     printf(" ...%s\n", str);
-  }else{
-    if(last_level == 1) printf("\n");
+  }
+  else {
+    if (last_level == 1) printf("\n");
     int level2 = level + (level_change == -1 ? -1 : 0);
-    for(int i = 0; i < level2; i++) printf("|  ");
+    for (int i = 0; i < level2; i++) printf("|  ");
     printf("%s%s", (level_change == -1 ? "..." : ""), str);
-    if(level_change != 1) printf("\n");
+    if (level_change != 1) printf("\n");
   }
   level     += level_change;
   last_level = level_change;
@@ -145,32 +168,31 @@ inline void messagePrintf(int level_change, char const* fmt, ...){
 /*!
  * @brief 將兩個字串用人類習慣的檔名排序方式排序
  *
- * 例如 a1 \< a2 \< a3 \< a10 \< a12 \< a20, 
+ * 例如 a1 \< a2 \< a3 \< a10 \< a12 \< a20,
  * 而不是 a1 \< a10 \< a12 \< a2 \< a20 \< a3
  *
  * @param [in] f1 第一個字串
  * @param [in] f2 第二個字串
  * @return \c  true/false 表 \b f1是否該排在f2前面
  */
-inline bool filenameCompare(std::string const& f1, std::string const& f2){
+inline bool filenameCompare(std::string const& f1, std::string const& f2) {
   char const* s1 = f1.c_str();
   char const* s2 = f2.c_str();
   int l1 = f1.length();
   int l2 = f2.length();
   int i1, i2;
-  for(i1 = i2 = 0; i1 < l1 || i2 < l2; i1++, i2++){
-    if(isdigit(s1[i1]) && isdigit(s2[i2])){
+  for (i1 = i2 = 0; i1 < l1 || i2 < l2; i1++, i2++) {
+    if (isdigit(s1[i1]) && isdigit(s2[i2])) {
       int n1 = atoi(s1 + i1);
       int n2 = atoi(s2 + i2);
-      if(n1 != n2){
+      if (n1 != n2)
         return (n1 < n2);
-      }
       while(i1 + 1 < l1 && isdigit(s1[i1 + 1])) i1++;
       while(i2 + 1 < l2 && isdigit(s2[i2 + 1])) i2++;
-    }else{
-      if(s1[i1] != s2[i2]){
+    }
+    else {
+      if(s1[i1] != s2[i2])
         return s1[i1] < s2[i2];
-      }
     }
   }
   return false;