From 7e95bada9aebfbfdcbdbdfcad350c5df520d86e1 Mon Sep 17 00:00:00 2001
From: martin <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 12:07:45 +0100
Subject: [PATCH 01/19] Add specific building instructions

---
 BUILD.md    | 15 +++++++++++++
 BUILD.md~   |  0
 README      |  3 +++
 README~     | 26 ++++++++++++++++++++++
 Readme.txt  | 14 +-----------
 Readme.txt~ | 62 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 6 files changed, 107 insertions(+), 13 deletions(-)
 create mode 100644 BUILD.md
 create mode 100644 BUILD.md~
 create mode 100644 README~
 create mode 100644 Readme.txt~

diff --git a/BUILD.md b/BUILD.md
new file mode 100644
index 00000000..da2bd8a8
--- /dev/null
+++ b/BUILD.md
@@ -0,0 +1,15 @@
+# Windows
+Windows => Use precompiled binary, or compile it with VS2008/2010 (Express or pro, Pro will allow you to enable Opemp in CMVS)
+        => Use CMake GUI in order to generate the Visual Studio project file (in ./program you will find the main CMakeLists.txt).
+
+# Linux compilations (Ubuntu used as example)
+```
+#Prepare and empty machine for building:
+sudo apt-get update -qq && sudo apt-get install -qq
+sudo apt-get -y install git jpeg boost boost-graph
+git clone https://github.com/pmoulon/CMVS-PMVS
+mkdir CMVS-PMVS_build && cd CMVS-PMVS_build
+cmake ../CMVS-PMVS/program
+make
+sudo make install
+```
diff --git a/BUILD.md~ b/BUILD.md~
new file mode 100644
index 00000000..e69de29b
diff --git a/README b/README
index 25a7d3ab..a192e1f5 100644
--- a/README
+++ b/README
@@ -1,5 +1,8 @@
 This is a modified version of CMVS/PMVS.
 
+Building:
+ - See [BUILD.md](https://github.com/pmoulon/CMVS-PMVS/BUILD.md)
+
 Main modification:
  - cross platform compilation Linux, Windows => CMake build system generator.
  - added bug fix from Nghia Ho.
diff --git a/README~ b/README~
new file mode 100644
index 00000000..295aef36
--- /dev/null
+++ b/README~
@@ -0,0 +1,26 @@
+This is a modified version of CMVS/PMVS.
+
+Building:
+ - See [BUILD.md](https://github.com/pmoulon/CMVS-PMVS)
+
+Main modification:
+ - cross platform compilation Linux, Windows => CMake build system generator.
+ - added bug fix from Nghia Ho.
+ - make PLY, PSET, PATCH export faster and optional.
+ - Replaced GSL simplex/lmfit with nlopt optimizer
+ - Replaced image loading routines with CImg. Now PPMs are supported properly, with optional support for PNG and TIFF
+ - Replaced BLAS/LAPACK with Eigen
+ - Updated internal jpeg library and miniBoost
+ - CMake-system now supports system boost, jpeg and other libraries if available. 
+ - Replaced pthread with tinycthread to get rid of pthread.dll on Windows
+
+Authors : 
+[Original code author]  Yasutaka Furukawa http://www.cs.washington.edu/homes/furukawa/
+
+[Initial Cmake multiplatform port ]	Pierre moulon pmoulon[AT]gmail.com
+
+--------------------
+- Web ressources : - 
+--------------------
+[CMVS/PMVS] http://http://grail.cs.washington.edu/software/cmvs/
+[CMake version] http://opensourcephotogrammetry.blogspot.com/ https://github.com/pmoulon/CMVS-PMVS
diff --git a/Readme.txt b/Readme.txt
index 6eeb2d45..96f87633 100644
--- a/Readme.txt
+++ b/Readme.txt
@@ -18,19 +18,7 @@ Date : 13 July 2011
 -- Compilation  --
 --------------------
 
-Windows => Use precompiled binary, or compile it with VS2008/2010 (Express or pro, Pro will allow you to enable Opemp in CMVS)
-        => Use CMake GUI in order to generate the Visual Studio project file (in ./program you will find the main CMakeLists.txt).
-
-Linux => use makefile in program/main.
-
- Or use CMake build system :
-=> Install the following libraries : jpeg boost boost-graph
- $ mkdir OutputLinux
- $ cd OutputLinux
- $ cmake . ..
- $ make
- => That's all. Openmp is not activated yet. Add openmp in the cmvs link option and define the _OPENMP cxx flags
-
+- See https://github.com/pmoulon/CMVS-PMVS/BUILD.md
 
 --------------------
 ----  Notes :   ----
diff --git a/Readme.txt~ b/Readme.txt~
new file mode 100644
index 00000000..6eeb2d45
--- /dev/null
+++ b/Readme.txt~
@@ -0,0 +1,62 @@
+Authors : 
+[main author] 		Yasutaka Furukawa furukawa[AT]cs.washington.edu
+[main author] 		Jean Ponce Jean.Ponce[AT]ens.fr
+[Windows Porting] 	Pierre moulon pmoulon[AT]gmail.com
+
+Special thanks to ASTRE Henri for the PTHREAD 64 bits lib and dll : http://www.visual-experiments.com/
+
+Date : 13 July 2011
+
+--------------------
+- Web ressources : - 
+--------------------
+[CMVS] http://grail.cs.washington.edu/software/cmvs/
+[windows porting] http://francemapping.free.fr/Portfolio/Prog3D/CMVS.html
+
+
+--------------------
+-- Compilation  --
+--------------------
+
+Windows => Use precompiled binary, or compile it with VS2008/2010 (Express or pro, Pro will allow you to enable Opemp in CMVS)
+        => Use CMake GUI in order to generate the Visual Studio project file (in ./program you will find the main CMakeLists.txt).
+
+Linux => use makefile in program/main.
+
+ Or use CMake build system :
+=> Install the following libraries : jpeg boost boost-graph
+ $ mkdir OutputLinux
+ $ cd OutputLinux
+ $ cmake . ..
+ $ make
+ => That's all. Openmp is not activated yet. Add openmp in the cmvs link option and define the _OPENMP cxx flags
+
+
+--------------------
+----  Notes :   ----
+--------------------
+To test it : ./WinX-VS2010/Readme.txt "Usage of the binaries"
+
+What have been done on native Yasutaka Furukawa source code :
+
+- Create the CMake build system for CMVS/PMVS2 and required library.
+
+- Optimize a little bit JPEG image loading.
+
+- Update CMVS source code in order to compile.
+
+- Add changes from Nghia Ho http://nghiaho.com/
+  - memoize pow(2,X)
+  - Change GSL simplex to lmfit.
+
+- Replaced GSL simplex/lmfit with nlopt optimizer
+
+- Replaced image loading routines with CImg. Now PPMs are supported properly, with optional support for PNG and TIFF
+
+- Replaced BLAS/LAPACK with Eigen
+
+- Updated internal jpeg library and miniBoost
+
+- CMake-system now supports system boost, jpeg and other libraries if available. 
+
+- Replaced pthread with tinycthread to get rid of pthread.dll on Windows

From 47ac49e929f2eb8db205db324bc0a38c5008ba27 Mon Sep 17 00:00:00 2001
From: martin <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 12:08:39 +0100
Subject: [PATCH 02/19] Add specific building instructions

---
 BUILD.md~   |  0
 README~     | 26 ----------------------
 Readme.txt~ | 62 -----------------------------------------------------
 3 files changed, 88 deletions(-)
 delete mode 100644 BUILD.md~
 delete mode 100644 README~
 delete mode 100644 Readme.txt~

diff --git a/BUILD.md~ b/BUILD.md~
deleted file mode 100644
index e69de29b..00000000
diff --git a/README~ b/README~
deleted file mode 100644
index 295aef36..00000000
--- a/README~
+++ /dev/null
@@ -1,26 +0,0 @@
-This is a modified version of CMVS/PMVS.
-
-Building:
- - See [BUILD.md](https://github.com/pmoulon/CMVS-PMVS)
-
-Main modification:
- - cross platform compilation Linux, Windows => CMake build system generator.
- - added bug fix from Nghia Ho.
- - make PLY, PSET, PATCH export faster and optional.
- - Replaced GSL simplex/lmfit with nlopt optimizer
- - Replaced image loading routines with CImg. Now PPMs are supported properly, with optional support for PNG and TIFF
- - Replaced BLAS/LAPACK with Eigen
- - Updated internal jpeg library and miniBoost
- - CMake-system now supports system boost, jpeg and other libraries if available. 
- - Replaced pthread with tinycthread to get rid of pthread.dll on Windows
-
-Authors : 
-[Original code author]  Yasutaka Furukawa http://www.cs.washington.edu/homes/furukawa/
-
-[Initial Cmake multiplatform port ]	Pierre moulon pmoulon[AT]gmail.com
-
---------------------
-- Web ressources : - 
---------------------
-[CMVS/PMVS] http://http://grail.cs.washington.edu/software/cmvs/
-[CMake version] http://opensourcephotogrammetry.blogspot.com/ https://github.com/pmoulon/CMVS-PMVS
diff --git a/Readme.txt~ b/Readme.txt~
deleted file mode 100644
index 6eeb2d45..00000000
--- a/Readme.txt~
+++ /dev/null
@@ -1,62 +0,0 @@
-Authors : 
-[main author] 		Yasutaka Furukawa furukawa[AT]cs.washington.edu
-[main author] 		Jean Ponce Jean.Ponce[AT]ens.fr
-[Windows Porting] 	Pierre moulon pmoulon[AT]gmail.com
-
-Special thanks to ASTRE Henri for the PTHREAD 64 bits lib and dll : http://www.visual-experiments.com/
-
-Date : 13 July 2011
-
---------------------
-- Web ressources : - 
---------------------
-[CMVS] http://grail.cs.washington.edu/software/cmvs/
-[windows porting] http://francemapping.free.fr/Portfolio/Prog3D/CMVS.html
-
-
---------------------
--- Compilation  --
---------------------
-
-Windows => Use precompiled binary, or compile it with VS2008/2010 (Express or pro, Pro will allow you to enable Opemp in CMVS)
-        => Use CMake GUI in order to generate the Visual Studio project file (in ./program you will find the main CMakeLists.txt).
-
-Linux => use makefile in program/main.
-
- Or use CMake build system :
-=> Install the following libraries : jpeg boost boost-graph
- $ mkdir OutputLinux
- $ cd OutputLinux
- $ cmake . ..
- $ make
- => That's all. Openmp is not activated yet. Add openmp in the cmvs link option and define the _OPENMP cxx flags
-
-
---------------------
-----  Notes :   ----
---------------------
-To test it : ./WinX-VS2010/Readme.txt "Usage of the binaries"
-
-What have been done on native Yasutaka Furukawa source code :
-
-- Create the CMake build system for CMVS/PMVS2 and required library.
-
-- Optimize a little bit JPEG image loading.
-
-- Update CMVS source code in order to compile.
-
-- Add changes from Nghia Ho http://nghiaho.com/
-  - memoize pow(2,X)
-  - Change GSL simplex to lmfit.
-
-- Replaced GSL simplex/lmfit with nlopt optimizer
-
-- Replaced image loading routines with CImg. Now PPMs are supported properly, with optional support for PNG and TIFF
-
-- Replaced BLAS/LAPACK with Eigen
-
-- Updated internal jpeg library and miniBoost
-
-- CMake-system now supports system boost, jpeg and other libraries if available. 
-
-- Replaced pthread with tinycthread to get rid of pthread.dll on Windows

From 79107f57bc7151fc7996fd2d4bd4deeb71a89f57 Mon Sep 17 00:00:00 2001
From: martin <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 12:10:10 +0100
Subject: [PATCH 03/19] Add specific building instructions

---
 README  |  2 +-
 README~ | 26 ++++++++++++++++++++++++++
 2 files changed, 27 insertions(+), 1 deletion(-)
 create mode 100644 README~

diff --git a/README b/README
index a192e1f5..9027231e 100644
--- a/README
+++ b/README
@@ -1,7 +1,7 @@
 This is a modified version of CMVS/PMVS.
 
 Building:
- - See [BUILD.md](https://github.com/pmoulon/CMVS-PMVS/BUILD.md)
+ - See https://github.com/pmoulon/CMVS-PMVS/BUILD.md
 
 Main modification:
  - cross platform compilation Linux, Windows => CMake build system generator.
diff --git a/README~ b/README~
new file mode 100644
index 00000000..a192e1f5
--- /dev/null
+++ b/README~
@@ -0,0 +1,26 @@
+This is a modified version of CMVS/PMVS.
+
+Building:
+ - See [BUILD.md](https://github.com/pmoulon/CMVS-PMVS/BUILD.md)
+
+Main modification:
+ - cross platform compilation Linux, Windows => CMake build system generator.
+ - added bug fix from Nghia Ho.
+ - make PLY, PSET, PATCH export faster and optional.
+ - Replaced GSL simplex/lmfit with nlopt optimizer
+ - Replaced image loading routines with CImg. Now PPMs are supported properly, with optional support for PNG and TIFF
+ - Replaced BLAS/LAPACK with Eigen
+ - Updated internal jpeg library and miniBoost
+ - CMake-system now supports system boost, jpeg and other libraries if available. 
+ - Replaced pthread with tinycthread to get rid of pthread.dll on Windows
+
+Authors : 
+[Original code author]  Yasutaka Furukawa http://www.cs.washington.edu/homes/furukawa/
+
+[Initial Cmake multiplatform port ]	Pierre moulon pmoulon[AT]gmail.com
+
+--------------------
+- Web ressources : - 
+--------------------
+[CMVS/PMVS] http://http://grail.cs.washington.edu/software/cmvs/
+[CMake version] http://opensourcephotogrammetry.blogspot.com/ https://github.com/pmoulon/CMVS-PMVS

From 97aeff11f96cda73f25e99661910af2468756639 Mon Sep 17 00:00:00 2001
From: martin <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 12:12:39 +0100
Subject: [PATCH 04/19] remove README~

---
 README~ | 26 --------------------------
 1 file changed, 26 deletions(-)
 delete mode 100644 README~

diff --git a/README~ b/README~
deleted file mode 100644
index a192e1f5..00000000
--- a/README~
+++ /dev/null
@@ -1,26 +0,0 @@
-This is a modified version of CMVS/PMVS.
-
-Building:
- - See [BUILD.md](https://github.com/pmoulon/CMVS-PMVS/BUILD.md)
-
-Main modification:
- - cross platform compilation Linux, Windows => CMake build system generator.
- - added bug fix from Nghia Ho.
- - make PLY, PSET, PATCH export faster and optional.
- - Replaced GSL simplex/lmfit with nlopt optimizer
- - Replaced image loading routines with CImg. Now PPMs are supported properly, with optional support for PNG and TIFF
- - Replaced BLAS/LAPACK with Eigen
- - Updated internal jpeg library and miniBoost
- - CMake-system now supports system boost, jpeg and other libraries if available. 
- - Replaced pthread with tinycthread to get rid of pthread.dll on Windows
-
-Authors : 
-[Original code author]  Yasutaka Furukawa http://www.cs.washington.edu/homes/furukawa/
-
-[Initial Cmake multiplatform port ]	Pierre moulon pmoulon[AT]gmail.com
-
---------------------
-- Web ressources : - 
---------------------
-[CMVS/PMVS] http://http://grail.cs.washington.edu/software/cmvs/
-[CMake version] http://opensourcephotogrammetry.blogspot.com/ https://github.com/pmoulon/CMVS-PMVS

From 60bb7eb836c34cbb422afa8eadc12ea92e928198 Mon Sep 17 00:00:00 2001
From: martin <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 12:14:01 +0100
Subject: [PATCH 05/19] Add openmp line

---
 BUILD.md | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/BUILD.md b/BUILD.md
index da2bd8a8..da862daf 100644
--- a/BUILD.md
+++ b/BUILD.md
@@ -12,4 +12,6 @@ mkdir CMVS-PMVS_build && cd CMVS-PMVS_build
 cmake ../CMVS-PMVS/program
 make
 sudo make install
+
+ => Openmp is not activated yet. Add openmp in the cmvs link option and define the _OPENMP cxx flags
 ```

From d8a0ac7d8928d46b0b2775486a321b1dea66ec4f Mon Sep 17 00:00:00 2001
From: mad-de <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 14:46:21 +0100
Subject: [PATCH 06/19] Update pmvs2.cc

---
 program/main/pmvs2.cc | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/program/main/pmvs2.cc b/program/main/pmvs2.cc
index 5aa4b42f..49cda9d9 100644
--- a/program/main/pmvs2.cc
+++ b/program/main/pmvs2.cc
@@ -8,7 +8,7 @@ using namespace std;
 
 int main(int argc, char* argv[]) {
   if (argc < 3) {
-    cerr << "Usage: " << argv[0] << " prefix option_file [Optional export]" << endl
+    cout << "Usage: " << argv[0] << " prefix option_file [Optional export]" << endl
          << endl
          << "--------------------------------------------------" << endl
          << "level       1    csize    2" << endl

From 4e0913f17f26372d1a5951b8d1b90dea2b90dc55 Mon Sep 17 00:00:00 2001
From: mad-de <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 14:50:10 +0100
Subject: [PATCH 07/19] Update option.cc

---
 program/base/pmvs/option.cc | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/program/base/pmvs/option.cc b/program/base/pmvs/option.cc
index 8c9d1d80..19a8d51f 100644
--- a/program/base/pmvs/option.cc
+++ b/program/base/pmvs/option.cc
@@ -117,28 +117,28 @@ void Soption::init(const std::string prefix, const std::string option) {
   if (m_useBound)
     initBindexes(sbimages);
 
-  cerr << "--------------------------------------------------" << endl  
+  cout << "--------------------------------------------------" << endl  
        << "--- Summary of specified options ---" << endl;
-  cerr << "# of timages: " << (int)m_timages.size();
+  cout << "# of timages: " << (int)m_timages.size();
   if (m_tflag == -1)
-    cerr << " (range specification)" << endl;
+    cout << " (range specification)" << endl;
   else
-    cerr << " (enumeration)" << endl;
-  cerr << "# of oimages: " << (int)m_oimages.size();
+    cout << " (enumeration)" << endl;
+  cout << "# of oimages: " << (int)m_oimages.size();
   if (m_oflag == -1)
-    cerr << " (range specification)" << endl;
+    cout << " (range specification)" << endl;
   else if (0 <= m_oflag)
-    cerr << " (enumeration)" << endl;
+    cout << " (enumeration)" << endl;
   else if (m_oflag == -2)
-    cerr << " (vis.dat is used)" << endl;
+    cout << " (vis.dat is used)" << endl;
   else if (m_oflag == -3)
-    cerr << " (not used)" << endl;
+    cout << " (not used)" << endl;
 
-  cerr << "level: " << m_level << "  csize: " << m_csize << endl
+  cout << "level: " << m_level << "  csize: " << m_csize << endl
        << "threshold: " << m_threshold << "  wsize: " << m_wsize << endl
        << "minImageNum: " << m_minImageNum << "  CPU: " << m_CPU << endl
        << "useVisData: " << m_useVisData << "  sequence: " << m_sequence << endl;
-  cerr << "--------------------------------------------------" << endl;
+  cout << "--------------------------------------------------" << endl;
 }
 
 void Soption::initOimages(void) {

From 19d46074a573c299579d905065cbe31287bdcf87 Mon Sep 17 00:00:00 2001
From: mad-de <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 14:59:14 +0100
Subject: [PATCH 08/19] Update photoSetS.cc

---
 program/base/image/photoSetS.cc | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/program/base/image/photoSetS.cc b/program/base/image/photoSetS.cc
index 9fb6717c..368e92b2 100644
--- a/program/base/image/photoSetS.cc
+++ b/program/base/image/photoSetS.cc
@@ -25,7 +25,7 @@ void CphotoSetS::init(const std::vector<int>& images, const std::string prefix,
   m_prefix = prefix;
   m_maxLevel = max(1, maxLevel);
   m_photos.resize(m_num);
-  cerr << "Reading images: " << flush;
+  cout << "Reading images: " << flush;
   for (int index = 0; index < m_num; ++index) {
     const int image = m_images[index];
 
@@ -56,7 +56,7 @@ void CphotoSetS::init(const std::vector<int>& images, const std::string prefix,
         m_photos[index].alloc();
       else
         m_photos[index].alloc(1);
-      cerr << '*' << flush;
+      cout << '*' << flush;
     }
     // try 4 digits
     else {
@@ -73,7 +73,7 @@ void CphotoSetS::init(const std::vector<int>& images, const std::string prefix,
         m_photos[index].alloc();
       else
         m_photos[index].alloc(1);
-      cerr << '*' << flush;
+      cout << '*' << flush;
     }
 
     /*
@@ -91,10 +91,10 @@ void CphotoSetS::init(const std::vector<int>& images, const std::string prefix,
       m_photos[index].alloc();
     else
       m_photos[index].alloc(1);
-    cerr << '*' << flush;
+    cout << '*' << flush;
     */
   }
-  cerr << endl;
+  cout << endl;
   const int margin = size / 2;
   m_size = 2 * margin + 1;
 }

From 4acf8002b93fdf81501b775fe02bb21483124583 Mon Sep 17 00:00:00 2001
From: mad-de <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 15:01:55 +0100
Subject: [PATCH 09/19] Update harris.cc

---
 program/base/pmvs/harris.cc | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/program/base/pmvs/harris.cc b/program/base/pmvs/harris.cc
index 4d68cc9a..ebd77637 100644
--- a/program/base/pmvs/harris.cc
+++ b/program/base/pmvs/harris.cc
@@ -166,7 +166,7 @@ void Charris::run(const std::vector<unsigned char>& image,
 		  const int gspeedup, const float sigma,
 		  std::multiset<Cpoint> & result) {
 
-  cerr << "Harris running ..." << flush;
+  cout << "Harris running ..." << flush;
   m_width = width;       m_height = height;
   m_sigmaD = sigma;      m_sigmaI = sigma;
   init(image, mask, edge);
@@ -219,5 +219,5 @@ void Charris::run(const std::vector<unsigned char>& image,
       }
     }
 
-  cerr << (int)result.size() << " harris done" << endl;
+  cout << (int)result.size() << " harris done" << endl;
 }

From cb40c86a7dbc771b95eaa583d453dbebfa44cabd Mon Sep 17 00:00:00 2001
From: mad-de <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 15:02:26 +0100
Subject: [PATCH 10/19] Update dog.cc

---
 program/base/pmvs/dog.cc | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/program/base/pmvs/dog.cc b/program/base/pmvs/dog.cc
index cd133668..040d43c8 100644
--- a/program/base/pmvs/dog.cc
+++ b/program/base/pmvs/dog.cc
@@ -119,7 +119,7 @@ void Cdog::run(const std::vector<unsigned char>& image,
 	       const float firstScale,   // 1.4f
 	       const float lastScale,    // 4.0f
 	       std::multiset<Cpoint> & result) {
-  cerr << "DoG running..." << flush;
+  cout << "DoG running..." << flush;
   m_width = width;
   m_height = height;
   
@@ -221,7 +221,7 @@ void Cdog::run(const std::vector<unsigned char>& image,
       }
     }
   
-  cerr << (int)result.size() << " dog done" << endl;  
+  cout << (int)result.size() << " dog done" << endl;  
 }
 
 void Cdog::setRes(const float sigma,

From bcbb4f5b2a9364b54b1aa0d4711d164a2df9d6cd Mon Sep 17 00:00:00 2001
From: martin <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 15:31:53 +0100
Subject: [PATCH 11/19] change cerr to cout where not necessary

---
 program/base/cmvs/bundle.cc          | 86 ++++++++++++++--------------
 program/base/pmvs/detectFeatures.cc  |  4 +-
 program/base/pmvs/expand.cc          |  8 +--
 program/base/pmvs/filter.cc          | 28 ++++-----
 program/base/pmvs/findMatch.cc       |  2 +-
 program/base/pmvs/option.cc          |  2 +-
 program/base/pmvs/patchOrganizerS.cc |  4 +-
 program/base/pmvs/seed.cc            | 12 ++--
 program/main/cmvs.cc                 |  2 +-
 program/main/genOption.cc            |  2 +-
 10 files changed, 75 insertions(+), 75 deletions(-)

diff --git a/program/base/cmvs/bundle.cc b/program/base/cmvs/bundle.cc
index 82f8c719..ceb2fc34 100644
--- a/program/base/cmvs/bundle.cc
+++ b/program/base/cmvs/bundle.cc
@@ -40,7 +40,7 @@ void Cbundle::prep(const std::string prefix, const int imageThreshold,
                    const float coverageThreshold,
                    const int pnumThreshold, const int CPU) {
   if (pnumThreshold != 0) {
-    cerr << "Should use pnumThreshold = 0" << endl;
+    cout << "Should use pnumThreshold = 0" << endl;
     exit (1);
   }
   
@@ -59,9 +59,9 @@ void Cbundle::prep(const std::string prefix, const int imageThreshold,
   
   char buffer[1024];
   sprintf(buffer, "%sbundle.rd.out", prefix.c_str());
-  cerr << "Reading bundle..." << flush;
+  cout << "Reading bundle..." << flush;
   readBundle(buffer);
-  cerr << endl;
+  cout << endl;
 
   vector<int> images;
   for (int c = 0; c < m_cnum; ++c)
@@ -71,9 +71,9 @@ void Cbundle::prep(const std::string prefix, const int imageThreshold,
   m_maxLevel = 12;
   m_pss.init(images, prefix, m_maxLevel + 1, 5, 0);
   
-  cerr << "Set widths/heights..." << flush;
+  cout << "Set widths/heights..." << flush;
   setWidthsHeightsLevels();
-  cerr << "done" << flush;
+  cout << "done" << flush;
 }
 
 void Cbundle::prep2(void) {
@@ -83,40 +83,40 @@ void Cbundle::prep2(void) {
     m_sfms2.resize((int)m_coords.size());
     startTimer();
     setScoreThresholds();
-    cerr << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
+    cout << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
     startTimer();
-    cerr << "slimNeighborsSetLinks..." << flush;
+    cout << "slimNeighborsSetLinks..." << flush;
     slimNeighborsSetLinks();
-    cerr << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
+    cout << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
   }
   
   // Improve visibility by using texture analysis
   startTimer();  
-  cerr << "mergeSFM..." << flush;
+  cout << "mergeSFM..." << flush;
   mergeSfMP();
-  cerr << '\t' << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
+  cout << '\t' << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
 
   m_sfms2.clear();
   m_sfms2.resize((int)m_coords.size());
-  cerr << "setScoreThresholds..." << flush;
+  cout << "setScoreThresholds..." << flush;
   startTimer();
   setScoreThresholds();
-  cerr << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
+  cout << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
 
   // Remove redundant images first
-  cerr << "sRemoveImages... " << flush;
+  cout << "sRemoveImages... " << flush;
   startTimer();
 
   sRemoveImages();
-  cerr << '\t' << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
+  cout << '\t' << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
 
   // use m_removed to change m_visibles and update m_neighbors
   startTimer();
   resetVisibles();
   setNeighbors();
-  cerr << "slimNeighborsSetLinks..." << flush;
+  cout << "slimNeighborsSetLinks..." << flush;
   slimNeighborsSetLinks();
-  cerr << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
+  cout << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
   
   // Init m_timages by mutually exclusive clustering
   setTimages();
@@ -132,12 +132,12 @@ void Cbundle::run(const std::string prefix, const int imageThreshold,
   prep(prefix, imageThreshold, tau, scoreRatioThreshold,
        coverageThreshold, pnumThreshold, CPU);
 
-  cerr << '\t' << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
+  cout << '\t' << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
 
   prep2();
   
   // Assumed variables that must be set properly here
-  cerr << "Adding images: " << endl;
+  cout << "Adding images: " << endl;
   startTimer();
   // Add images
   // Repeat until all the clusters become at most m_imageThreshold.
@@ -170,7 +170,7 @@ void Cbundle::run(const std::string prefix, const int imageThreshold,
     if (change == 0)
       break;
   }
-  cerr << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
+  cout << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
   
   m_oimages.resize((int)m_timages.size());
 
@@ -265,28 +265,28 @@ void Cbundle::sRemoveImages(void) {
   const int tenth = max(1, m_cnum / 10);
   for (int i = 0; i < (int)vvi.size(); ++i) {
     if (i % tenth == 0)
-      cerr << '*' << flush;
+      cout << '*' << flush;
     const int image = vvi[i][1];
     checkImage(image);
   }
-  cerr << endl;
+  cout << endl;
     
-  cerr << "Kept: ";
+  cout << "Kept: ";
   int kept = 0;
   for (int c = 0; c < m_cnum; ++c)
     if (m_removed[c] == 0) {
       ++kept;
-      cerr << c << ' ';
+      cout << c << ' ';
     }
-  cerr << endl << endl;
+  cout << endl << endl;
   
-  cerr << "Removed: ";
+  cout << "Removed: ";
   for (int c = 0; c < m_cnum; ++c)
     if (m_removed[c]) {
-      cerr << c << ' ';
+      cout << c << ' ';
     }
-  cerr << endl;
-  cerr << "sRemoveImages: " << m_cnum << " -> " << kept << flush;
+  cout << endl;
+  cout << "sRemoveImages: " << m_cnum << " -> " << kept << flush;
 }
 
 void Cbundle::resetVisibles(void) {
@@ -448,10 +448,10 @@ void Cbundle::setTimages(void) {
   else
     divideImages(lhs, m_timages);
   
-  cerr << endl << "Cluster sizes: " << endl;
+  cout << endl << "Cluster sizes: " << endl;
   for (int i = 0; i < (int)m_timages.size(); ++i)
-    cerr << (int)m_timages[i].size() << ' ';
-  cerr << endl;
+    cout << (int)m_timages[i].size() << ' ';
+  cout << endl;
 }
 
 void Cbundle::divideImages(const std::vector<int>& lhs,
@@ -558,7 +558,7 @@ void Cbundle::readBundle(const std::string file) {
   int cnum, pnum;
   ifstr >> cnum >> pnum;
   vector<int> ids;        ids.resize(cnum);
-  cerr << cnum << " cameras -- " << pnum << " points in bundle file" << endl;
+  cout << cnum << " cameras -- " << pnum << " points in bundle file" << endl;
   m_cnum = 0;
   for (int c = 0; c < cnum; ++c) {
     ids[c] = -1;
@@ -577,7 +577,7 @@ void Cbundle::readBundle(const std::string file) {
   const int tenth = max(1, pnum / 10);
   for (int p = 0; p < pnum; ++p) {
     if (p % tenth == 0)
-      cerr << '*' << flush;
+      cout << '*' << flush;
     int num;    Vec3f color;
     Vec4f coord;
     ifstr >> coord[0] >> coord[1] >> coord[2]
@@ -621,7 +621,7 @@ void Cbundle::readBundle(const std::string file) {
   ifstr.close();
   setNeighbors();
   
-  cerr << endl << m_cnum << " cameras -- " << m_pnum << " points" << flush;
+  cout << endl << m_cnum << " cameras -- " << m_pnum << " points" << flush;
 }
 
 void Cbundle::findPNeighbors(sfcnn<const float*, 3, float>& tree,
@@ -673,7 +673,7 @@ void Cbundle::mergeSfMPThread(void) {
     if (pid != -1 && m_merged[pid])
       pid = -2;
     if (m_count % tenth == 0)
-      cerr << '*' << flush;
+      cout << '*' << flush;
     ++m_count;
     mtx_unlock(&m_lock);
     if (pid == -2)
@@ -786,9 +786,9 @@ void Cbundle::mergeSfMP(void) {
   }
   
   // Based on m_puf, reset m_coords, m_coords, m_visibles, m_vpoints
-  cerr << "resetPoints..." << flush;
+  cout << "resetPoints..." << flush;
   resetPoints();
-  cerr << "done" << endl;
+  cout << "done" << endl;
 
   delete m_puf;
   m_puf = NULL;
@@ -796,7 +796,7 @@ void Cbundle::mergeSfMP(void) {
   m_ptree = NULL;
    
   const int npnum = (int)m_coords.size();
-  cerr << "Rep counts: " << cpnum << " -> " << npnum << "  " << flush;
+  cout << "Rep counts: " << cpnum << " -> " << npnum << "  " << flush;
 }
 
 // Based on m_puf, reset m_coords, m_coords, m_visibles, m_vpoints
@@ -907,7 +907,7 @@ void Cbundle::setCluster(const int p) {
     }
     // If none of the visibles images are
     if (find == 0) {
-      cerr << "Impossible in setscoresclustersthread" << endl
+      cout << "Impossible in setscoresclustersthread" << endl
            << (int)m_visibles[p].size() << endl
            << (int)m_sfms2[p].m_satisfied << endl;
       exit (1);
@@ -998,8 +998,8 @@ void Cbundle::addImagesP(void) {
   }
 
   for (int i = 0; i < (int)m_addnums.size(); ++i)
-    cerr << m_addnums[i] << ' ';
-  cerr << endl;
+    cout << m_addnums[i] << ' ';
+  cout << endl;
   
   int totalnum = 0;
   for (int c = 0; c < (int)m_timages.size(); ++c)
@@ -1009,7 +1009,7 @@ void Cbundle::addImagesP(void) {
     if (m_removed[c] == 0)
       ++beforenum;
   
-  cerr << "Image nums: "
+  cout << "Image nums: "
        << m_cnum << " -> " << beforenum <<  " -> " << totalnum << endl;
 }
 
@@ -1389,7 +1389,7 @@ void Cbundle::writeVis(void) {
   ofstr.close();
 
 
-  cerr << numer / (float)denom << " images in vis on the average" << endl;
+  cout << numer / (float)denom << " images in vis on the average" << endl;
   
 }
 
diff --git a/program/base/pmvs/detectFeatures.cc b/program/base/pmvs/detectFeatures.cc
index 5c7012f4..abb86b40 100644
--- a/program/base/pmvs/detectFeatures.cc
+++ b/program/base/pmvs/detectFeatures.cc
@@ -39,7 +39,7 @@ void CdetectFeatures::run(const CphotoSetS& pss, const int num,
   for (int i = 0; i < m_CPU; ++i)
     thrd_join(threads[i], NULL);
   //----------------------------------------------------------------------
-  cerr << "done" << endl;
+  cout << "done" << endl;
 }
 
 int CdetectFeatures::runThreadTmp(void* arg) {
@@ -61,7 +61,7 @@ void CdetectFeatures::runThread(void) {
       break;
     
     const int image = m_ppss->m_images[index];
-    cerr << image << ' ' << flush;
+    cout << image << ' ' << flush;
 
     //?????????????  May need file lock, because targetting images
     //should not overlap among multiple processors.    
diff --git a/program/base/pmvs/expand.cc b/program/base/pmvs/expand.cc
index 76d32554..3e86b325 100644
--- a/program/base/pmvs/expand.cc
+++ b/program/base/pmvs/expand.cc
@@ -43,24 +43,24 @@ void Cexpand::run(void) {
   // set queue
   m_fm.m_pos.collectPatches(m_queue);
 
-  cerr << "Expanding patches..." << flush;
+  cout << "Expanding patches..." << flush;
   vector<thrd_t> threads(m_fm.m_CPU);
   for (int c = 0; c < m_fm.m_CPU; ++c)
     thrd_create(&threads[c], &expandThreadTmp, (void*)this);
   for (int c = 0; c < m_fm.m_CPU; ++c)
     thrd_join(threads[c], NULL); 
   
-  cerr << endl
+  cout << endl
        << "---- EXPANSION: " << (time(NULL) - starttime) << " secs ----" << endl;
 
   const int trial = accumulate(m_ecounts.begin(), m_ecounts.end(), 0);
   const int fail0 = accumulate(m_fcounts0.begin(), m_fcounts0.end(), 0);
   const int fail1 = accumulate(m_fcounts1.begin(), m_fcounts1.end(), 0);
   const int pass = accumulate(m_pcounts.begin(), m_pcounts.end(), 0);
-  cerr << "Total pass fail0 fail1 refinepatch: "
+  cout << "Total pass fail0 fail1 refinepatch: "
        << trial << ' ' << pass << ' '
        << fail0 << ' ' << fail1 << ' ' << pass + fail1 << endl;
-  cerr << "Total pass fail0 fail1 refinepatch: "
+  cout << "Total pass fail0 fail1 refinepatch: "
        << 100 * trial / (float)trial << ' '
        << 100 * pass / (float)trial << ' '
        << 100 * fail0 / (float)trial << ' '
diff --git a/program/base/pmvs/filter.cc b/program/base/pmvs/filter.cc
index d2faeeb2..e4dcd6fd 100644
--- a/program/base/pmvs/filter.cc
+++ b/program/base/pmvs/filter.cc
@@ -36,14 +36,14 @@ void Cfilter::filterOutside(void) {
   time_t tv;
   time(&tv); 
   time_t curtime = tv;
-  cerr << "FilterOutside" << endl;
+  cout << "FilterOutside" << endl;
   //??? notice (1) here to avoid removing m_fix=1
   m_fm.m_pos.collectPatches(1);
 
   const int psize = (int)m_fm.m_pos.m_ppatches.size();  
   m_gains.resize(psize);
   
-  cerr << "mainbody: " << flush;
+  cout << "mainbody: " << flush;
   
   m_fm.m_count = 0;
   vector<thrd_t> threads(m_fm.m_CPU);
@@ -51,7 +51,7 @@ void Cfilter::filterOutside(void) {
     thrd_create(&threads[i], &filterOutsideThreadTmp, (void*)this);
   for (int i = 0; i < m_fm.m_CPU; ++i)
     thrd_join(threads[i], NULL);
-  cerr << endl;
+  cout << endl;
 
   // delete patches with positive m_gains
   int count = 0;
@@ -76,10 +76,10 @@ void Cfilter::filterOutside(void) {
   ave /= denom;
   ave2 /= denom;
   ave2 = sqrt(max(0.0, ave2 - ave * ave));
-  cerr << "Gain (ave/var): " << ave << ' ' << ave2 << endl;
+  cout << "Gain (ave/var): " << ave << ' ' << ave2 << endl;
   
   time(&tv);
-  cerr << (int)m_fm.m_pos.m_ppatches.size() << " -> "
+  cout << (int)m_fm.m_pos.m_ppatches.size() << " -> "
        << (int)m_fm.m_pos.m_ppatches.size() - count << " ("
        << 100 * ((int)m_fm.m_pos.m_ppatches.size() - count) / (float)m_fm.m_pos.m_ppatches.size()
        << "%)\t" << (tv - curtime) / CLOCKS_PER_SEC << " secs" << endl;
@@ -216,7 +216,7 @@ void Cfilter::filterExact(void) {
   time_t tv;
   time(&tv); 
   time_t curtime = tv;
-  cerr << "Filter Exact: " << flush;
+  cout << "Filter Exact: " << flush;
 
   //??? cannot use (1) because we use patch.m_id to set newimages,....
   m_fm.m_pos.collectPatches();
@@ -234,7 +234,7 @@ void Cfilter::filterExact(void) {
     thrd_create(&threads0[i], &filterExactThreadTmp, (void*)this);  
   for (int i = 0; i < m_fm.m_CPU; ++i)
     thrd_join(threads0[i], NULL);
-  cerr << endl;
+  cout << endl;
 
   //----------------------------------------------------------------------
   for (int p = 0; p < psize; ++p) {
@@ -291,7 +291,7 @@ void Cfilter::filterExact(void) {
     }
   }
   time(&tv); 
-  cerr << (int)m_fm.m_pos.m_ppatches.size() << " -> "
+  cout << (int)m_fm.m_pos.m_ppatches.size() << " -> "
        << (int)m_fm.m_pos.m_ppatches.size() - count << " ("
        << 100 * ((int)m_fm.m_pos.m_ppatches.size() - count) / (float)m_fm.m_pos.m_ppatches.size()
        << "%)\t" << (tv - curtime) / CLOCKS_PER_SEC << " secs" << endl;
@@ -312,7 +312,7 @@ void Cfilter::filterExactThread(void) {
     if (m_fm.m_tnum <= image)
       break;
     
-    cerr << '*' << flush;
+    cout << '*' << flush;
     
     const int& w = m_fm.m_pos.m_gwidths[image];
     const int& h = m_fm.m_pos.m_gheights[image];
@@ -516,7 +516,7 @@ void Cfilter::filterNeighbor(const int times) {
   time_t tv;
   time(&tv); 
   time_t curtime = tv;
-  cerr << "FilterNeighbor:\t" << flush;
+  cout << "FilterNeighbor:\t" << flush;
 
   //??? notice (1) to avoid removing m_fix=1
   m_fm.m_pos.collectPatches(1);
@@ -558,7 +558,7 @@ void Cfilter::filterNeighbor(const int times) {
     }
   }
   time(&tv);
-  cerr << (int)m_fm.m_pos.m_ppatches.size() << " -> "
+  cout << (int)m_fm.m_pos.m_ppatches.size() << " -> "
        << (int)m_fm.m_pos.m_ppatches.size() - count << " ("
        << 100 * ((int)m_fm.m_pos.m_ppatches.size() - count) / (float)m_fm.m_pos.m_ppatches.size()
        << "%)\t" << (tv - curtime) / CLOCKS_PER_SEC << " secs" << endl;
@@ -571,7 +571,7 @@ void Cfilter::filterSmallGroups(void) {
   time_t tv;
   time(&tv); 
   time_t curtime = tv;
-  cerr << "FilterGroups:\t" << flush;
+  cout << "FilterGroups:\t" << flush;
   m_fm.m_pos.collectPatches();
   if (m_fm.m_pos.m_ppatches.empty())
     return;
@@ -619,7 +619,7 @@ void Cfilter::filterSmallGroups(void) {
   }
   
   const int threshold = max(20, psize / 10000);
-  cerr << threshold << endl;
+  cout << threshold << endl;
   
   bite = size.begin();
   eite = size.end();
@@ -651,7 +651,7 @@ void Cfilter::filterSmallGroups(void) {
     ++bpatch;
   }
   time(&tv);
-  cerr << (int)m_fm.m_pos.m_ppatches.size() << " -> "
+  cout << (int)m_fm.m_pos.m_ppatches.size() << " -> "
        << (int)m_fm.m_pos.m_ppatches.size() - count << " ("
        << 100 * ((int)m_fm.m_pos.m_ppatches.size() - count) / (float)m_fm.m_pos.m_ppatches.size()
        << "%)\t" << (tv - curtime)/CLOCKS_PER_SEC << " secs" << endl;
diff --git a/program/base/pmvs/findMatch.cc b/program/base/pmvs/findMatch.cc
index 5f434432..2e803998 100644
--- a/program/base/pmvs/findMatch.cc
+++ b/program/base/pmvs/findMatch.cc
@@ -225,7 +225,7 @@ void CfindMatch::run(void) {
     ++m_depth;
   }
   time(&tv);
-  cerr << "---- Total: " << (tv - curtime)/CLOCKS_PER_SEC << " secs ----" << endl;
+  cout << "---- Total: " << (tv - curtime)/CLOCKS_PER_SEC << " secs ----" << endl;
 }
 
 void CfindMatch::write(const std::string prefix, bool bExportPLY, bool bExportPatch, bool bExportPSet) {
diff --git a/program/base/pmvs/option.cc b/program/base/pmvs/option.cc
index 19a8d51f..fb6f5001 100644
--- a/program/base/pmvs/option.cc
+++ b/program/base/pmvs/option.cc
@@ -284,7 +284,7 @@ void Soption::initBindexes(const std::string sbimages) {
     exit (1);
   }
   
-  cerr << "Reading bimages" << endl;
+  cout << "Reading bimages" << endl;
   int itmp;
   ifstr >> itmp;
   for (int i = 0; i < itmp; ++i) {
diff --git a/program/base/pmvs/patchOrganizerS.cc b/program/base/pmvs/patchOrganizerS.cc
index dd11bcf6..0e82c613 100644
--- a/program/base/pmvs/patchOrganizerS.cc
+++ b/program/base/pmvs/patchOrganizerS.cc
@@ -132,7 +132,7 @@ void CpatchOrganizerS::readPatches(void) {
 
     string header;    int pnum;
     ifstr >> header >> pnum;
-    cerr << image << ' ' << pnum << " patches" << endl;
+    cout << image << ' ' << pnum << " patches" << endl;
     for (int p = 0; p < pnum; ++p) {
       Ppatch ppatch(new Cpatch());
       ifstr >> *ppatch;
@@ -171,7 +171,7 @@ void CpatchOrganizerS::readPatches(void) {
     
     string header;    int pnum;
     ifstr >> header >> pnum;
-    cerr << image << ' ' << pnum << " patches" << endl;
+    cout << image << ' ' << pnum << " patches" << endl;
     
     for (int p = 0; p < pnum; ++p) {
       Ppatch ppatch(new Cpatch());
diff --git a/program/base/pmvs/seed.cc b/program/base/pmvs/seed.cc
index 65f827aa..cc78963a 100644
--- a/program/base/pmvs/seed.cc
+++ b/program/base/pmvs/seed.cc
@@ -58,7 +58,7 @@ void Cseed::run(void) {
   random_shuffle(vitmp.begin(), vitmp.end());
   m_fm.m_jobs.insert(m_fm.m_jobs.end(), vitmp.begin(), vitmp.end());
 
-  cerr << "adding seeds " << endl;
+  cout << "adding seeds " << endl;
   
   m_fm.m_pos.clearCounts();
 
@@ -79,18 +79,18 @@ void Cseed::run(void) {
   for (int i = 0; i < m_fm.m_CPU; ++i)
     thrd_join(threads[i], NULL);
   //----------------------------------------------------------------------
-  cerr << "done" << endl;
+  cout << "done" << endl;
   time(&tv);
-  cerr << "---- Initial: " << (tv - curtime)/CLOCKS_PER_SEC << " secs ----" << endl;
+  cout << "---- Initial: " << (tv - curtime)/CLOCKS_PER_SEC << " secs ----" << endl;
 
   const int trial = accumulate(m_scounts.begin(), m_scounts.end(), 0);
   const int fail0 = accumulate(m_fcounts0.begin(), m_fcounts0.end(), 0);
   const int fail1 = accumulate(m_fcounts1.begin(), m_fcounts1.end(), 0);
   const int pass = accumulate(m_pcounts.begin(), m_pcounts.end(), 0);
-  cerr << "Total pass fail0 fail1 refinepatch: "
+  cout << "Total pass fail0 fail1 refinepatch: "
        << trial << ' ' << pass << ' '
        << fail0 << ' ' << fail1 << ' ' << pass + fail1 << endl;
-  cerr << "Total pass fail0 fail1 refinepatch: "
+  cout << "Total pass fail0 fail1 refinepatch: "
        << 100 * trial / (float)trial << ' '
        << 100 * pass / (float)trial << ' '
        << 100 * fail0 / (float)trial << ' '
@@ -198,7 +198,7 @@ void Cseed::initialMatch(const int index, const int id) {
       }
     }
   }
-  cerr << '(' << index << ',' << totalcount << ')' << flush;
+  cout << '(' << index << ',' << totalcount << ')' << flush;
 }
 
 void Cseed::collectCells(const int index0, const int index1,
diff --git a/program/main/cmvs.cc b/program/main/cmvs.cc
index ca996ac1..22605bd9 100644
--- a/program/main/cmvs.cc
+++ b/program/main/cmvs.cc
@@ -8,7 +8,7 @@ using namespace std;
 
 int main(int argc, char* argv[]) {
   if (argc < 2) {
-    cerr << "Usage: " << argv[0] << " prefix maximage[=100] CPU[=4]" << endl
+    cout << "Usage: " << argv[0] << " prefix maximage[=100] CPU[=4]" << endl
          << endl
          << "You should choose maximage based on the amount of memory in your machine." << endl
          << "CPU should be the number of (virtual) CPUs or cores in your machine." << endl
diff --git a/program/main/genOption.cc b/program/main/genOption.cc
index 6ba24457..b0c4abe0 100644
--- a/program/main/genOption.cc
+++ b/program/main/genOption.cc
@@ -9,7 +9,7 @@ using namespace std;
 
 int main(int argc, char* argv[]) {
   if (argc < 2) {
-    cerr << "Usage: " << argv[0]
+    cout << "Usage: " << argv[0]
          << " prefix level[=1] csize[=2] threshold[=0.7] wsize[=7]"
          << " minImageNum[=3] CPU[=8]" << endl
          << endl

From 2abf85b8c25ec0a069dfeb853396ab630ccf6160 Mon Sep 17 00:00:00 2001
From: martin <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 15:39:50 +0100
Subject: [PATCH 12/19] prepare to pr

---
 README     |  3 +++
 Readme.txt | 14 +++++++++++++-
 2 files changed, 16 insertions(+), 1 deletion(-)

diff --git a/README b/README
index 9027231e..b8a7e26d 100644
--- a/README
+++ b/README
@@ -1,8 +1,11 @@
 This is a modified version of CMVS/PMVS.
 
+<<<<<<< HEAD
 Building:
  - See https://github.com/pmoulon/CMVS-PMVS/BUILD.md
 
+=======
+>>>>>>> parent of 7e95bad... Add specific building instructions
 Main modification:
  - cross platform compilation Linux, Windows => CMake build system generator.
  - added bug fix from Nghia Ho.
diff --git a/Readme.txt b/Readme.txt
index 96f87633..6eeb2d45 100644
--- a/Readme.txt
+++ b/Readme.txt
@@ -18,7 +18,19 @@ Date : 13 July 2011
 -- Compilation  --
 --------------------
 
-- See https://github.com/pmoulon/CMVS-PMVS/BUILD.md
+Windows => Use precompiled binary, or compile it with VS2008/2010 (Express or pro, Pro will allow you to enable Opemp in CMVS)
+        => Use CMake GUI in order to generate the Visual Studio project file (in ./program you will find the main CMakeLists.txt).
+
+Linux => use makefile in program/main.
+
+ Or use CMake build system :
+=> Install the following libraries : jpeg boost boost-graph
+ $ mkdir OutputLinux
+ $ cd OutputLinux
+ $ cmake . ..
+ $ make
+ => That's all. Openmp is not activated yet. Add openmp in the cmvs link option and define the _OPENMP cxx flags
+
 
 --------------------
 ----  Notes :   ----

From 1de37cc379490c5f7506c2e959bc4d05e5103294 Mon Sep 17 00:00:00 2001
From: martin <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 15:40:18 +0100
Subject: [PATCH 13/19] Revert "Add specific building instructions"

This reverts commit 79107f57bc7151fc7996fd2d4bd4deeb71a89f57.
---
 README | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README b/README
index b8a7e26d..68d9cde6 100644
--- a/README
+++ b/README
@@ -2,7 +2,7 @@ This is a modified version of CMVS/PMVS.
 
 <<<<<<< HEAD
 Building:
- - See https://github.com/pmoulon/CMVS-PMVS/BUILD.md
+ - See [BUILD.md](https://github.com/pmoulon/CMVS-PMVS/BUILD.md)
 
 =======
 >>>>>>> parent of 7e95bad... Add specific building instructions

From 6f779a558ab08655d0e3b083740bd52378723e49 Mon Sep 17 00:00:00 2001
From: martin <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 15:46:07 +0100
Subject: [PATCH 14/19] prepare to pr

---
 BUILD.md | 17 -----------------
 README   |  6 ------
 2 files changed, 23 deletions(-)
 delete mode 100644 BUILD.md

diff --git a/BUILD.md b/BUILD.md
deleted file mode 100644
index da862daf..00000000
--- a/BUILD.md
+++ /dev/null
@@ -1,17 +0,0 @@
-# Windows
-Windows => Use precompiled binary, or compile it with VS2008/2010 (Express or pro, Pro will allow you to enable Opemp in CMVS)
-        => Use CMake GUI in order to generate the Visual Studio project file (in ./program you will find the main CMakeLists.txt).
-
-# Linux compilations (Ubuntu used as example)
-```
-#Prepare and empty machine for building:
-sudo apt-get update -qq && sudo apt-get install -qq
-sudo apt-get -y install git jpeg boost boost-graph
-git clone https://github.com/pmoulon/CMVS-PMVS
-mkdir CMVS-PMVS_build && cd CMVS-PMVS_build
-cmake ../CMVS-PMVS/program
-make
-sudo make install
-
- => Openmp is not activated yet. Add openmp in the cmvs link option and define the _OPENMP cxx flags
-```
diff --git a/README b/README
index 68d9cde6..25a7d3ab 100644
--- a/README
+++ b/README
@@ -1,11 +1,5 @@
 This is a modified version of CMVS/PMVS.
 
-<<<<<<< HEAD
-Building:
- - See [BUILD.md](https://github.com/pmoulon/CMVS-PMVS/BUILD.md)
-
-=======
->>>>>>> parent of 7e95bad... Add specific building instructions
 Main modification:
  - cross platform compilation Linux, Windows => CMake build system generator.
  - added bug fix from Nghia Ho.

From 7b7b5c42c0fe303dfb43ff64eb698430494f6f6e Mon Sep 17 00:00:00 2001
From: martin <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 15:57:43 +0100
Subject: [PATCH 15/19] revert one cerr

---
 program/base/cmvs/bundle.cc | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/program/base/cmvs/bundle.cc b/program/base/cmvs/bundle.cc
index ceb2fc34..a687e735 100644
--- a/program/base/cmvs/bundle.cc
+++ b/program/base/cmvs/bundle.cc
@@ -907,7 +907,7 @@ void Cbundle::setCluster(const int p) {
     }
     // If none of the visibles images are
     if (find == 0) {
-      cout << "Impossible in setscoresclustersthread" << endl
+      cerr << "Impossible in setscoresclustersthread" << endl
            << (int)m_visibles[p].size() << endl
            << (int)m_sfms2[p].m_satisfied << endl;
       exit (1);

From 40c011224f2bd4927d6f7db8d94f4d5e4241bce4 Mon Sep 17 00:00:00 2001
From: martin <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 15:59:31 +0100
Subject: [PATCH 16/19] revert one cerr

---
 program/base/cmvs/bundle.cc | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/program/base/cmvs/bundle.cc b/program/base/cmvs/bundle.cc
index a687e735..d0be6d08 100644
--- a/program/base/cmvs/bundle.cc
+++ b/program/base/cmvs/bundle.cc
@@ -40,7 +40,7 @@ void Cbundle::prep(const std::string prefix, const int imageThreshold,
                    const float coverageThreshold,
                    const int pnumThreshold, const int CPU) {
   if (pnumThreshold != 0) {
-    cout << "Should use pnumThreshold = 0" << endl;
+    cerr << "Should use pnumThreshold = 0" << endl;
     exit (1);
   }
   

From 76e8bff6fe1eeb335b2f0646399d956c8732c19d Mon Sep 17 00:00:00 2001
From: martin <cev_madde@msn.com>
Date: Tue, 16 Feb 2016 18:43:55 +0100
Subject: [PATCH 17/19] Add Openmp to CMaklists.txt

---
 program/CMakeLists.txt | 27 +++++++++++++++++++++++++++
 1 file changed, 27 insertions(+)

diff --git a/program/CMakeLists.txt b/program/CMakeLists.txt
index 2a4cd941..61d2d438 100644
--- a/program/CMakeLists.txt
+++ b/program/CMakeLists.txt
@@ -30,6 +30,9 @@ INCLUDE_DIRECTORIES(
 FIND_PACKAGE(Threads REQUIRED)
 SET(PMVS_LIBRARIES ${PMVS_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT})
 
+# OPTIONS
+OPTION(CIMG_USE_OPENMP "Enable OpenMP parallelization in cimg" ON)
+
 # Eigen
 SET( EIGEN3_INCLUDE_DIR ${CMAKE_SOURCE_DIR}/thirdParty/Eigen )
 SET( CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${EIGEN3_INCLUDE_DIR}/cmake" )
@@ -95,6 +98,30 @@ IF(PMVS_USE_TIFF)
   ENDIF(TIFF_FOUND)
 ENDIF(PMVS_USE_TIFF)
 
+# OpenMP
+# ==============================================================================
+# OpenMP detection
+# ==============================================================================
+IF(CIMG_USE_OPENMP)
+  FIND_PACKAGE(OpenMP)
+  IF(OPENMP_FOUND)
+    SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+    OPTION(CIMG_USE_OPENMP "Use OpenMP for parallelization in cimg" ON)
+    ADD_DEFINITIONS(-DCIMG_USE_OPENMP)
+    IF (NOT MSVC)
+      IF ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
+        # for those using the clang with OpenMP support
+        LIST(APPEND OPENMVG_LIBRARY_DEPENDENCIES iomp)
+      ELSE()
+        LIST(APPEND OPENMVG_LIBRARY_DEPENDENCIES gomp)
+      ENDIF()
+    ENDIF(NOT MSVC)
+  ENDIF(OPENMP_FOUND)
+ELSE(CIMG_USE_OPENMP)
+    OPTION(OpenMVG_USE_OPENMP "Use OpenMP for parallelization in cimg" OFF)
+    UPDATE_CACHE_VARIABLE(CIMG_USE_OPENMP OFF)
+    REMOVE_DEFINITIONS(-DCIMG_USE_OPENMP)
+ENDIF(CIMG_USE_OPENMP)
 
 ADD_SUBDIRECTORY(base)
 ADD_SUBDIRECTORY(main)

From 6ddad19c68d1e0e8094628748af4c6282c95076c Mon Sep 17 00:00:00 2001
From: martin <cev_madde@msn.com>
Date: Thu, 18 Feb 2016 00:30:05 +0100
Subject: [PATCH 18/19] cerr breaks at this point

---
 0                            |    0
 program/base/cmvs/bundle.cc  |    2 +-
 program/base/cmvs/bundle.cc~ | 1476 ++++++++++++++++++++++++++++++++++
 3 files changed, 1477 insertions(+), 1 deletion(-)
 create mode 100644 0
 create mode 100644 program/base/cmvs/bundle.cc~

diff --git a/0 b/0
new file mode 100644
index 00000000..e69de29b
diff --git a/program/base/cmvs/bundle.cc b/program/base/cmvs/bundle.cc
index d0be6d08..cf11a82d 100644
--- a/program/base/cmvs/bundle.cc
+++ b/program/base/cmvs/bundle.cc
@@ -515,7 +515,7 @@ void Cbundle::divideImages(const std::vector<int>& lhs,
     }
     
     if (cand1.empty() || cand2.empty()) {
-      cerr << "Error. Normalized cuts produced an empty cluster: "
+      cout << "Error. Normalized cuts produced an empty cluster: "
            << (int)part.size() << " -> "
            << (int)cand1.size() << ' '
            << (int)cand2.size() << endl;
diff --git a/program/base/cmvs/bundle.cc~ b/program/base/cmvs/bundle.cc~
new file mode 100644
index 00000000..d0be6d08
--- /dev/null
+++ b/program/base/cmvs/bundle.cc~
@@ -0,0 +1,1476 @@
+#include <fstream>
+#include <iterator>
+#include <numeric> //PM
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+#include "graclus.h"
+#include "bundle.h"
+#define _USE_MATH_DEFINES
+#include <math.h>
+
+extern "C"
+{
+int boundary_points;
+int spectral_initialization = 0;
+int cutType;
+int memory_saving;
+};
+
+using namespace std;
+using namespace boost;
+using namespace CMVS;
+
+Cbundle::Cbundle(void) {
+  m_CPU = 8;
+  m_junit = 100;
+  mtx_init(&m_lock, mtx_plain | mtx_recursive);
+  m_debug = 0;
+  m_puf = NULL;
+  m_ptree = NULL;
+}
+
+Cbundle::~Cbundle() {
+}
+
+void Cbundle::prep(const std::string prefix, const int imageThreshold,
+                   const int tau, const float scoreRatioThreshold,
+                   const float coverageThreshold,
+                   const int pnumThreshold, const int CPU) {
+  if (pnumThreshold != 0) {
+    cerr << "Should use pnumThreshold = 0" << endl;
+    exit (1);
+  }
+  
+  m_prefix = prefix;
+  m_imageThreshold = imageThreshold;
+  m_tau = tau;
+  m_scoreRatioThreshold = scoreRatioThreshold;
+  m_coverageThreshold = coverageThreshold;
+  m_pnumThreshold = pnumThreshold;
+  m_CPU = CPU;
+  
+  m_linkThreshold = 2.0f;
+
+  m_dscale = 1 / 100.0f;
+  m_dscale2 = 1.0f;
+  
+  char buffer[1024];
+  sprintf(buffer, "%sbundle.rd.out", prefix.c_str());
+  cout << "Reading bundle..." << flush;
+  readBundle(buffer);
+  cout << endl;
+
+  vector<int> images;
+  for (int c = 0; c < m_cnum; ++c)
+    images.push_back(c);
+
+  m_dlevel = 7;
+  m_maxLevel = 12;
+  m_pss.init(images, prefix, m_maxLevel + 1, 5, 0);
+  
+  cout << "Set widths/heights..." << flush;
+  setWidthsHeightsLevels();
+  cout << "done" << flush;
+}
+
+void Cbundle::prep2(void) {
+  // Used in mergeSfMP now.
+  {
+    m_pweights.resize((int)m_coords.size(), 1);
+    m_sfms2.resize((int)m_coords.size());
+    startTimer();
+    setScoreThresholds();
+    cout << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
+    startTimer();
+    cout << "slimNeighborsSetLinks..." << flush;
+    slimNeighborsSetLinks();
+    cout << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
+  }
+  
+  // Improve visibility by using texture analysis
+  startTimer();  
+  cout << "mergeSFM..." << flush;
+  mergeSfMP();
+  cout << '\t' << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
+
+  m_sfms2.clear();
+  m_sfms2.resize((int)m_coords.size());
+  cout << "setScoreThresholds..." << flush;
+  startTimer();
+  setScoreThresholds();
+  cout << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
+
+  // Remove redundant images first
+  cout << "sRemoveImages... " << flush;
+  startTimer();
+
+  sRemoveImages();
+  cout << '\t' << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
+
+  // use m_removed to change m_visibles and update m_neighbors
+  startTimer();
+  resetVisibles();
+  setNeighbors();
+  cout << "slimNeighborsSetLinks..." << flush;
+  slimNeighborsSetLinks();
+  cout << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
+  
+  // Init m_timages by mutually exclusive clustering
+  setTimages();
+  m_oimages.resize((int)m_timages.size());  
+}
+
+void Cbundle::run(const std::string prefix, const int imageThreshold,
+                  const int tau, const float scoreRatioThreshold,
+                  const float coverageThreshold,
+                  const int pnumThreshold, const int CPU) {
+  startTimer();
+  
+  prep(prefix, imageThreshold, tau, scoreRatioThreshold,
+       coverageThreshold, pnumThreshold, CPU);
+
+  cout << '\t' << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
+
+  prep2();
+  
+  // Assumed variables that must be set properly here
+  cout << "Adding images: " << endl;
+  startTimer();
+  // Add images
+  // Repeat until all the clusters become at most m_imageThreshold.
+  while (1) {
+    addImagesP();
+    
+    int change = 0;
+    vector<vector<int> > newtimages;
+    cout << "Divide: " << flush;
+    for (int i = 0; i < (int)m_timages.size(); ++i) {
+      if ((int)m_timages[i].size() <= m_imageThreshold) {        
+        newtimages.push_back(m_timages[i]);
+        continue;
+      }
+      else {
+        cout << i << ' ';
+        
+        change = 1;
+        // divide
+        vector<vector<int> > vvi;
+        divideImages(m_timages[i], vvi);
+        for (int j = 0; j < (int)vvi.size(); ++j)
+          newtimages.push_back(vvi[j]);
+      }
+    }
+
+    cout << endl;
+    
+    m_timages.swap(newtimages);
+    if (change == 0)
+      break;
+  }
+  cout << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
+  
+  m_oimages.resize((int)m_timages.size());
+
+  writeCameraCenters();
+  
+  // Output results
+  writeVis();
+  writeGroups();
+}
+
+float Cbundle::computeLink(const int image0, const int image1) {
+  vector<int> common;
+  set_intersection(m_vpoints[image0].begin(), m_vpoints[image0].end(),
+                   m_vpoints[image1].begin(), m_vpoints[image1].end(),
+                   back_inserter(common));
+
+  float score = 0.0f;
+  for (int i = 0; i < (int)common.size(); ++i) {
+    const int pid = common[i];
+    vector<int> vtmp;
+    vtmp.push_back(image0);    vtmp.push_back(image1);
+    const float ftmp = computeScore2(m_coords[pid], vtmp);
+    if (m_sfms2[pid].m_score != 0.0f)
+      score += m_pweights[pid] *
+        ftmp / (m_sfms2[pid].m_scoreThreshold / m_scoreRatioThreshold);
+  }
+
+  return score;
+}
+
+void Cbundle::slimNeighborsSetLinks(void) {
+  const int maxneighbor = 30;
+  m_links.clear();
+  m_links.resize(m_cnum);
+
+#pragma omp parallel for
+  for (int image = 0; image < m_cnum; ++image) {
+    m_links[image].resize((int)m_neighbors[image].size());
+    
+    for (int i = 0; i < (int)m_neighbors[image].size(); ++i)
+      m_links[image][i] = computeLink(image, m_neighbors[image][i]);
+
+    if ((int)m_neighbors[image].size() < 2)
+      continue;
+
+    vector<int> newneighbors;
+    vector<float> newlinks;
+
+    vector<Vec2f> vv;
+    for (int i = 0; i < (int)m_neighbors[image].size(); ++i)
+      vv.push_back(Vec2(-m_links[image][i], m_neighbors[image][i]));
+    sort(vv.begin(), vv.end(), Svec2cmp<float>());
+
+    const int itmp = min(maxneighbor, (int)vv.size());
+    for (int i = 0; i < itmp; ++i) {
+      newneighbors.push_back((int)vv[i][1]);
+      newlinks.push_back(-vv[i][0]);
+    }
+    
+    m_neighbors[image].swap(newneighbors);
+    m_links[image].swap(newlinks);
+  }
+}
+
+void Cbundle::setScoreThresholds(void) {
+#pragma omp parallel for
+  for (int p = 0; p < (int)m_coords.size(); ++p) {
+    m_sfms2[p].m_scoreThreshold =
+      computeScore2(m_coords[p], m_visibles[p], m_sfms2[p].m_uimages)
+      * m_scoreRatioThreshold;
+  }
+}
+
+void Cbundle::sRemoveImages(void) {
+  m_removed.clear();
+  m_removed.resize(m_cnum, 0);
+    
+  m_allows.resize(m_cnum);
+  for (int c = 0; c < m_cnum; ++c)
+    m_allows[c] = (int)ceil((int)m_vpoints[c].size() *
+                            (1.0f - m_coverageThreshold));
+  
+  // Sort all the images in an increasing order of resolution
+  vector<Vec2i> vvi;
+  for (int c = 0; c < m_cnum; ++c) {
+    const int res =
+      m_pss.m_photos[c].getWidth(0) * m_pss.m_photos[c].getHeight(0);
+    vvi.push_back(Vec2i(res, c));
+  }
+  sort(vvi.begin(), vvi.end(), Svec2cmp<int>());
+  
+  const int tenth = max(1, m_cnum / 10);
+  for (int i = 0; i < (int)vvi.size(); ++i) {
+    if (i % tenth == 0)
+      cout << '*' << flush;
+    const int image = vvi[i][1];
+    checkImage(image);
+  }
+  cout << endl;
+    
+  cout << "Kept: ";
+  int kept = 0;
+  for (int c = 0; c < m_cnum; ++c)
+    if (m_removed[c] == 0) {
+      ++kept;
+      cout << c << ' ';
+    }
+  cout << endl << endl;
+  
+  cout << "Removed: ";
+  for (int c = 0; c < m_cnum; ++c)
+    if (m_removed[c]) {
+      cout << c << ' ';
+    }
+  cout << endl;
+  cout << "sRemoveImages: " << m_cnum << " -> " << kept << flush;
+}
+
+void Cbundle::resetVisibles(void) {
+  // reset m_visibles. remove "removed images" from the list.
+  for (int p = 0; p < (int)m_visibles.size(); ++p) {
+    vector<int> newimages;
+
+    setNewImages(p, -1, newimages);
+    m_visibles[p].swap(newimages);
+  }
+}
+
+void Cbundle::setNewImages(const int pid, const int rimage,
+                           std::vector<int>& newimages) {
+  for (int i = 0; i < (int)m_visibles[pid].size(); ++i) {
+    const int image = m_visibles[pid][i];
+    if (m_removed[image] || image == rimage)
+      continue;
+    newimages.push_back(image);
+  }
+}
+
+void Cbundle::checkImage(const int image) {
+  // For each SfM, check if it will be unsatisfied by removing image.
+  //  0: unsatisfy
+  //  1: satisfy->satisfy
+  //  2: satisfy->unsatisfy
+  m_statsT.resize((int)m_vpoints[image].size());
+  m_jobs.clear();
+  
+#pragma omp parallel for
+  for (int p = 0; p < (int)m_statsT.size(); ++p) {
+    const int pid = m_vpoints[image][p];
+    if (m_sfms2[pid].m_satisfied == 0)
+      m_statsT[p] = 0;
+    else {
+      m_statsT[p] = 1;
+
+      // if m_uimages of p is not removed, if image is not in
+      // m_uimages, the point is still satisfied.
+      int valid = 1;
+      int inside = 0;
+
+      for (int i = 0; i < (int)m_sfms2[pid].m_uimages.size(); ++i) {
+        const int itmp = m_sfms2[pid].m_uimages[i];
+        if (itmp == image)
+          inside = 1;
+        if (m_removed[itmp]) {
+          valid = 0;
+          break;
+        }
+      }
+      
+      if (valid == 1 && inside == 0)
+        continue;
+      
+      vector<int> newimages;
+      setNewImages(pid, image, newimages);
+      const float cscore = computeScore2(m_coords[pid], newimages);
+      if (cscore < m_sfms2[pid].m_scoreThreshold)
+        m_statsT[p] = 2;
+    }
+  }
+  
+  // For each image, how many SFM points are removed if you remove
+  // "image"
+  vector<int> decrements;
+  decrements.resize(m_cnum, 0);
+  
+  // Look at points with m_stastT[p] = 2, to see if m_allows are still
+  // above thresholds.
+  for (int p = 0; p < (int)m_statsT.size(); ++p) {
+    if (m_statsT[p] == 2) {
+      const int pid = m_vpoints[image][p];
+      for (int i = 0; i < (int)m_visibles[pid].size(); ++i) {
+        const int itmp = m_visibles[pid][i];
+        ++decrements[itmp];
+      }
+    }
+  }
+
+  // If m_allows can cover decrements, go ahead.
+  int rflag = 1;
+  for (int c = 0; c < m_cnum; ++c)
+    if (m_allows[c] < decrements[c]) {
+      rflag = 0;
+      break;
+    }
+
+  // remove
+  if (rflag) {
+    m_removed[image] = 1;
+    for (int p = 0; p < (int)m_statsT.size(); ++p)
+      if (m_statsT[p] == 2)
+        m_sfms2[m_vpoints[image][p]].m_satisfied = 0;
+
+    for (int c = 0; c < m_cnum; ++c)
+      m_allows[c] -= decrements[c];
+
+    // Update uimages for points that are still satisfied and still
+    // contains image in m_uimages
+#pragma omp parallel for
+    for (int p = 0; p < (int)m_statsT.size(); ++p) {
+      const int pid = m_vpoints[image][p];
+      if (m_statsT[p] == 1) {
+        int contain = 0;
+        for (int i = 0; i < (int)m_sfms2[pid].m_uimages.size(); ++i)
+          if (m_sfms2[pid].m_uimages[i] == image) {
+            contain = 1;
+            break;
+          }
+
+        if (contain) {
+          vector<int> newimages;
+          setNewImages(pid, -1, newimages);
+          const float cscore =
+            computeScore2(m_coords[pid], newimages, m_sfms2[pid].m_uimages);
+          if (cscore < m_sfms2[pid].m_scoreThreshold)
+            m_sfms2[pid].m_satisfied = 0;
+        }
+      }
+    }
+  }
+}
+
+void Cbundle::setNeighbors(void) {
+  m_neighbors.clear();
+  m_neighbors.resize(m_cnum);
+#pragma omp parallel for
+  for (int image = 0; image < m_cnum; ++image) {
+    vector<int> narray;
+    narray.resize(m_cnum, 0);
+    
+    for (int p = 0; p < (int)m_visibles.size(); ++p) {
+      if (!binary_search(m_visibles[p].begin(), m_visibles[p].end(),
+                         image))
+        continue;
+      
+      for (int i = 0; i < (int)m_visibles[p].size(); ++i)
+        narray[m_visibles[p][i]] = 1;
+    }
+    
+    for (int i = 0; i < m_cnum; ++i)
+      if (narray[i] && i != image)
+        m_neighbors[image].push_back(i);
+  } 
+}
+
+void Cbundle::setTimages(void) {
+  vector<int> lhs;
+  for (int c = 0; c < m_cnum; ++c)
+    if (m_removed[c] == 0)
+      lhs.push_back(c);
+
+  m_timages.clear();
+
+  if ((int)lhs.size() <= m_imageThreshold)
+    m_timages.push_back(lhs);
+  else
+    divideImages(lhs, m_timages);
+  
+  cout << endl << "Cluster sizes: " << endl;
+  for (int i = 0; i < (int)m_timages.size(); ++i)
+    cout << (int)m_timages[i].size() << ' ';
+  cout << endl;
+}
+
+void Cbundle::divideImages(const std::vector<int>& lhs,
+                           std::vector<std::vector<int> >& rhs) {
+  const float iratio = 125.0f / 150.0f;
+  // Candidates
+  list<vector<int> > candidates;
+  // initialize first cluster
+  candidates.push_back(vector<int>());
+  vector<int>& vitmp = candidates.back();
+  for (int i = 0; i < (int)lhs.size(); ++i)
+    vitmp.push_back(lhs[i]);
+  
+  // Iterate
+  while (1) {
+    if (candidates.empty())
+      break;
+    vector<int> cand = candidates.front();
+    candidates.pop_front();
+    if ((int)cand.size() <= m_imageThreshold * iratio) {
+      rhs.push_back(cand);
+      continue;
+    }
+    // Divide into 2 groups
+    vector<idxtype> xadj, adjncy, adjwgt, part;
+    const int nparts = 2;
+    const int cutType = 0;
+    // Define graphs
+    for (int i = 0; i < (int)cand.size(); ++i) {
+      xadj.push_back((int)adjncy.size());
+
+      for (int j = 0; j < (int)cand.size(); ++j) {
+        if (i == j)
+          continue;
+        // Check if cand[i] and cand[j] are connected
+        const int cid0 = cand[i];
+        const int cid1 = cand[j];
+        vector<int>::const_iterator nite = 
+          find(m_neighbors[cid0].begin(), m_neighbors[cid0].end(), cid1);
+        
+        if (nite != m_neighbors[cid0].end()) {
+          adjncy.push_back(j);
+          const int offset = nite - m_neighbors[cid0].begin();
+          const int weight =
+            min(5000, (int)floor(10.0f * m_links[cid0][offset] + 0.5f));
+          adjwgt.push_back(weight);
+        }
+      }
+    }
+    xadj.push_back((int)adjncy.size());
+
+    Cgraclus::runE(xadj, adjncy, adjwgt, nparts, cutType, part);
+    
+    // Divide into 2 groups
+    vector<int> cand1, cand2;
+    for (int i = 0; i < (int)part.size(); ++i) {
+      if (part[i] == 0)
+        cand1.push_back(cand[i]);
+      else
+        cand2.push_back(cand[i]);
+    }
+    
+    if (cand1.empty() || cand2.empty()) {
+      cerr << "Error. Normalized cuts produced an empty cluster: "
+           << (int)part.size() << " -> "
+           << (int)cand1.size() << ' '
+           << (int)cand2.size() << endl;
+      exit (1);
+    }
+    
+    if ((int)cand1.size() <= m_imageThreshold * iratio)
+      rhs.push_back(cand1);
+    else {
+      candidates.push_back(vector<int>());
+      (candidates.back()).swap(cand1);
+    }
+    if ((int)cand2.size() <= m_imageThreshold * iratio)
+      rhs.push_back(cand2);
+    else {
+      candidates.push_back(vector<int>());
+      (candidates.back()).swap(cand2);
+    }
+  }
+}
+
+void Cbundle::readBundle(const std::string file) {
+  // For each valid image, a list of connected images, and the second
+  // value is the number of common points.
+  ifstream ifstr;  ifstr.open(file.c_str());
+  if (!ifstr.is_open()) {
+    cerr << "Bundle file not found: " << file << endl;
+    exit (1);
+  }
+  while (1) {
+    unsigned char uctmp;
+    ifstr.read((char*)&uctmp, sizeof(unsigned char));
+    ifstr.putback(uctmp);
+    if (uctmp == '#') {
+      char buffer[1024];      ifstr.getline(buffer, 1024);
+    }
+    else
+      break;
+  }
+  int cnum, pnum;
+  ifstr >> cnum >> pnum;
+  vector<int> ids;        ids.resize(cnum);
+  cout << cnum << " cameras -- " << pnum << " points in bundle file" << endl;
+  m_cnum = 0;
+  for (int c = 0; c < cnum; ++c) {
+    ids[c] = -1;
+    float params[15];
+    for (int i = 0; i < 15; ++i)
+      ifstr >> params[i];
+    if (params[0] != 0.0f)
+      ids[c] = m_cnum++;
+  }
+
+  m_vpoints.resize(m_cnum);
+
+  m_coords.clear();  m_visibles.clear();
+  m_colors.clear();
+  m_pnum = 0;
+  const int tenth = max(1, pnum / 10);
+  for (int p = 0; p < pnum; ++p) {
+    if (p % tenth == 0)
+      cout << '*' << flush;
+    int num;    Vec3f color;
+    Vec4f coord;
+    ifstr >> coord[0] >> coord[1] >> coord[2]
+          >> color >> num;
+    coord[3] = 1.0f;
+    
+    vector<int> visibles;
+    for (int i = 0; i < num; ++i) {
+      int itmp;      ifstr >> itmp;
+      if (cnum <= itmp) {
+        double dtmp;
+        ifstr >> itmp >> dtmp >> dtmp;
+        continue;
+      }
+      if (ids[itmp] == -1) {
+        cerr << "impossible " << itmp << ' ' << ids[itmp] << endl;
+        exit (1);
+      }
+      visibles.push_back(ids[itmp]);
+
+      // Based on the bundler version, the number of parameters here
+      // are either 1 or 3. Currently, it seems to be 3.
+      double dtmp;
+      ifstr >> itmp >> dtmp >> dtmp;
+      //ifstr >> dtmp;
+    }
+    if ((int)visibles.size() < 2)
+      continue;
+
+    sort(visibles.begin(), visibles.end());
+    
+    for (int i = 0; i < (int)visibles.size(); ++i) {
+       m_vpoints[visibles[i]].push_back(m_pnum);
+    }
+    
+    m_visibles.push_back(visibles);
+    m_coords.push_back(coord);
+    m_colors.push_back(color);
+    ++m_pnum;
+  }
+  ifstr.close();
+  setNeighbors();
+  
+  cout << endl << m_cnum << " cameras -- " << m_pnum << " points" << flush;
+}
+
+void Cbundle::findPNeighbors(sfcnn<const float*, 3, float>& tree,
+                             const int pid, std::vector<int>& pneighbors) {
+  const float thresh = m_dscale2 * m_dscale2 *
+    m_minScales[pid] * m_minScales[pid];
+
+  vector<long unsigned int> ids;
+  vector<double> dists;
+  int k = min((int)m_coords.size(), 400);
+
+  while (1) {
+    ids.clear();
+    dists.clear();
+    tree.ksearch(&m_coords[pid][0], k, ids, dists);
+
+    if (thresh < dists[k - 1])
+      break;
+    if (k == (int)m_coords.size())
+      break;
+    
+    k = min((int)m_coords.size(), k * 2);
+  }
+
+  for (int i = 0; i < k; ++i) {
+    if (thresh < dists[i])
+      break;
+
+    // If distance from pid2 is also within threshold ok
+    const float thresh2 = m_dscale2 * m_dscale2 *
+      m_minScales[ids[i]] * m_minScales[ids[i]];
+    if (thresh2 < dists[i])
+      continue;
+    
+    if (ids[i] != (unsigned int)pid)
+      pneighbors.push_back(ids[i]);
+  }
+}
+
+void Cbundle::mergeSfMPThread(void) {
+  const int tenth = (int)m_coords.size() / 10;
+  while (1) {
+    int pid = -1;
+    mtx_lock(&m_lock);
+    if (!m_jobs.empty()) {
+      pid = m_jobs.front();
+      m_jobs.pop_front();
+    }
+    if (pid != -1 && m_merged[pid])
+      pid = -2;
+    if (m_count % tenth == 0)
+      cout << '*' << flush;
+    ++m_count;
+    mtx_unlock(&m_lock);
+    if (pid == -2)
+      continue;
+    if (pid == -1)
+      break;
+
+    // Process, and check later if merged flag is on
+    vector<int> pneighbors;
+    findPNeighbors(*m_ptree, pid, pneighbors);
+    const int psize = (int)pneighbors.size();
+    // visible images and its neighbors
+    vector<int> visibles = m_visibles[pid];
+    vector<int> vitmp;
+    for (int i = 0; i < (int)m_visibles[pid].size(); ++i) {
+      const int imagetmp = m_visibles[pid][i];
+      vitmp.clear();
+      mymerge(visibles, m_neighbors[imagetmp], vitmp);
+      vitmp.swap(visibles);
+    }
+    
+    vector<char> visflag(psize, 0);
+    // test visibility
+    for (int i = 0; i < psize; ++i) {
+      const int& pid2 = pneighbors[i];
+      if (my_isIntersect(visibles, m_visibles[pid2]))
+        visflag[i] = 1;
+    }
+    
+    // Now lock and try to register
+    mtx_lock(&m_lock);
+    // If the main one is removed, over... waste.
+    if (m_merged[pid] == 0) {
+      m_merged[pid] = 1;
+      for (int i = 0; i < psize; ++i) {
+        const int pid2 = pneighbors[i];
+        if (visflag[i] && m_merged[pid2] == 0) {
+          m_merged[pid2] = 1;
+          m_puf->union_set(pid, pid2);
+        }
+      }
+    }
+    mtx_unlock(&m_lock);
+  }
+}
+
+int Cbundle::mergeSfMPThreadTmp(void* arg) {
+  ((Cbundle*)arg)->mergeSfMPThread();
+  return 0;
+}
+
+void Cbundle::mergeSfMP(void) {
+  // Repeat certain times until no change
+  const int cpnum = (int)m_coords.size();
+  m_minScales.resize(cpnum);
+  for (int p = 0; p < cpnum; ++p) {
+    m_minScales[p] = INT_MAX/2;
+    for (int i = 0; i < (int)m_visibles[p].size(); ++i) {
+      const int image = m_visibles[p][i];
+      const float stmp =
+        m_pss.m_photos[image].getScale(m_coords[p], m_levels[image]);
+      m_minScales[p] = min(m_minScales[p], stmp);
+    }
+  }
+  
+  m_puf = new disjoint_sets_with_storage<>(cpnum);
+  for (int p = 0; p < cpnum; ++p)
+    m_puf->make_set(p);
+
+  {
+    // kdtree
+    vector<const float*> ppoints;
+    ppoints.resize((int)m_coords.size());
+    for (int i = 0; i < (int)m_coords.size(); ++i)
+      ppoints[i] = &(m_coords[i][0]);
+    m_ptree =
+      new sfcnn<const float*, 3, float>(&ppoints[0], (int)ppoints.size());
+    
+    m_merged.resize((int)m_coords.size(), 0);
+    m_jobs.clear();
+    vector<int> order;
+    order.resize(cpnum);
+    for (int p = 0; p < cpnum; ++p)
+      order[p] = p;
+    random_shuffle(order.begin(), order.end());
+    
+    for (int p = 0; p < cpnum; ++p)
+      m_jobs.push_back(order[p]);
+    
+    m_count = 0;
+    vector<thrd_t> threads(m_CPU);
+    for (int c = 0; c < m_CPU; ++c)
+      thrd_create(&threads[c], &mergeSfMPThreadTmp, (void*)this);
+    for (int c = 0; c < m_CPU; ++c)
+      thrd_join(threads[c], NULL);
+
+    int newpnum = 0;
+    // Mapping from m_pnum to new id for reps
+    vector<int> dict;
+    vector<int> reps;
+    dict.resize((int)m_coords.size(), -1);
+    for (int p = 0; p < (int)m_coords.size(); ++p) {
+      const int ptmp = m_puf->find_set(p);
+      if (p == ptmp) {
+        dict[p] = newpnum;
+        reps.push_back(p);
+        ++newpnum;
+      }
+    }
+  }
+  
+  // Based on m_puf, reset m_coords, m_coords, m_visibles, m_vpoints
+  cout << "resetPoints..." << flush;
+  resetPoints();
+  cout << "done" << endl;
+
+  delete m_puf;
+  m_puf = NULL;
+  delete m_ptree;
+  m_ptree = NULL;
+   
+  const int npnum = (int)m_coords.size();
+  cout << "Rep counts: " << cpnum << " -> " << npnum << "  " << flush;
+}
+
+// Based on m_puf, reset m_coords, m_coords, m_visibles, m_vpoints
+void Cbundle::resetPoints(void) {
+  vector<int> counts;
+  vector<int> smallestids;
+  counts.resize((int)m_coords.size(), 0);
+  smallestids.resize((int)m_coords.size(), INT_MAX/2);
+  for (int p = 0; p < (int)m_coords.size(); ++p) {
+    const int ptmp = m_puf->find_set(p);
+    smallestids[ptmp] = min(smallestids[ptmp], p);
+    ++counts[ptmp];
+  }
+  const int mthreshold = 2;
+  vector<Vec2i> vv2;
+  for (int p = 0; p < (int)m_coords.size(); ++p)
+    if (mthreshold <= counts[p])
+      vv2.push_back(Vec2i(smallestids[p], p));
+  sort(vv2.begin(), vv2.end(), Svec2cmp<int>());
+  
+  vector<int> newpids;
+  newpids.resize((int)m_coords.size(), -1);
+  int newpnum = (int)vv2.size();
+  for (int i = 0; i < newpnum; ++i)
+    newpids[vv2[i][1]] = i;
+  
+  vector<Vec4f> newcoords;  newcoords.resize(newpnum);
+  vector<vector<int> > newvisibles;  newvisibles.resize(newpnum);
+
+  for (int p = 0; p < (int)m_coords.size(); ++p) {
+    const int ptmp = m_puf->find_set(p);
+    if (counts[ptmp] < mthreshold)
+      continue;
+    
+    const int newpid = newpids[ptmp];
+    newcoords[newpid] += m_coords[p];
+    vector<int> vitmp;
+    mymerge(newvisibles[newpid], m_visibles[p], vitmp);
+    vitmp.swap(newvisibles[newpid]);
+  }  
+  
+  for (int i = 0; i < newpnum; ++i)
+    newcoords[i] = newcoords[i] / newcoords[i][3];
+  
+  // Update m_vpoints
+  m_coords.swap(newcoords);
+  m_visibles.swap(newvisibles);
+  
+  for (int c = 0; c < m_cnum; ++c)
+    m_vpoints[c].clear();  
+  
+  for (int p = 0; p < (int)m_coords.size(); ++p) {
+    for (int i = 0; i < (int)m_visibles[p].size(); ++i) {
+      const int itmp = m_visibles[p][i];
+      m_vpoints[itmp].push_back(p);
+    }
+  }
+
+  m_pweights.clear();
+  for (int i = 0; i < newpnum; ++i)
+    m_pweights.push_back(counts[vv2[i][1]]);
+}
+
+void Cbundle::setScoresClusters(void) {
+#pragma omp parallel for
+  for (int p = 0; p < (int)m_coords.size(); ++p) {
+    // if m_satisfied is 0, no hope (even all the images are in the
+    // cluster, not satisfied
+    if (m_sfms2[p].m_satisfied == 0)
+      continue;
+    
+    m_sfms2[p].m_satisfied = 2;
+    setCluster(p);
+  }
+}
+
+void Cbundle::setCluster(const int p) {
+  m_sfms2[p].m_score = -1.0f;
+  m_sfms2[p].m_cluster = -1;
+  for (int c = 0; c < (int)m_timages.size(); ++c) {
+    // Select images in cluster
+    vector<int> vitmp;
+    set_intersection(m_timages[c].begin(), m_timages[c].end(),
+                     m_visibles[p].begin(), m_visibles[p].end(),
+                     back_inserter(vitmp));
+    
+    const float stmp = computeScore2(m_coords[p], vitmp);
+    if (m_sfms2[p].m_score < stmp) {
+      m_sfms2[p].m_cluster = c;
+      m_sfms2[p].m_score = stmp;
+    }
+  }
+
+  // If no cluster contains 2 images
+  if (m_sfms2[p].m_cluster == -1) {
+    int find = 0;
+    for (int j = 0; j < (int)m_visibles[p].size(); ++j) {
+      if (find)
+        break;
+      for (int c = 0; c < (int)m_timages.size(); ++c)
+        if (binary_search(m_timages[c].begin(), m_timages[c].end(),
+                          m_visibles[p][j])) {
+          m_sfms2[p].m_cluster = c;
+          m_sfms2[p].m_score = 0.0f;
+          find = 1;
+          break;
+        }
+    }
+    // If none of the visibles images are
+    if (find == 0) {
+      cerr << "Impossible in setscoresclustersthread" << endl
+           << (int)m_visibles[p].size() << endl
+           << (int)m_sfms2[p].m_satisfied << endl;
+      exit (1);
+    }
+  }
+
+  if (m_sfms2[p].m_scoreThreshold <= m_sfms2[p].m_score) {
+    // SATISFIED
+    m_sfms2[p].m_satisfied = 1;
+
+    //  update m_lacks
+    mtx_lock(&m_lock);    
+    for (int i = 0; i < (int)m_visibles[p].size(); ++i) {
+      const int image = m_visibles[p][i];
+      --m_lacks[image];
+    }
+    mtx_unlock(&m_lock);
+  }
+}
+
+float Cbundle::angleScore(const Vec4f& ray0, const Vec4f& ray1) {
+  const static float lsigma = 5.0f * M_PI / 180.f;
+  const static float rsigma = 15.0f * M_PI / 180.0f;
+  const static float lsigma2 = 2.0f * lsigma * lsigma;
+  const static float rsigma2 = 2.0f * rsigma * rsigma;
+  const static float pivot = 20.0f * M_PI / 180.0f;
+
+  const float angle = acos(min(1.0f, ray0 * ray1));
+  const float diff = angle - pivot;
+
+  if (angle < pivot)
+    return exp(- diff * diff / lsigma2);
+  else
+    return exp(- diff * diff / rsigma2);
+}
+
+void Cbundle::setClusters(void) {
+#pragma omp parallel for
+  for (int p = 0; p < (int)m_sfms2.size(); ++p) {
+    if (m_sfms2[p].m_satisfied != 2)
+      continue;
+
+    setCluster(p);
+  }
+}
+
+void Cbundle::addImagesP(void) {
+  // set m_lacks (how many more sfm points need to be satisfied)
+  m_lacks.resize(m_cnum);
+  for (int c = 0; c < m_cnum; ++c) {
+    if (m_removed[c])
+      m_lacks[c] = 0;
+    else
+      m_lacks[c] =
+        (int)floor((int)m_vpoints[c].size() * m_coverageThreshold);
+  }
+
+  // Compute best score given all the images. Upper bound on the score
+  setScoresClusters();
+  
+  // Add images to clusters to make m_lacks at most 0 for all the
+  // images. In practice, for each cluster, identify corresponding sfm
+  // points that have not been satisfied. These points compute gains.
+  
+  // set m_adds for each sfm point, and add images
+
+  m_addnums.clear();
+  m_addnums.resize((int)m_timages.size(), 0);
+  
+  while (1) {
+    const int totalnum = addImages();
+    // Very end
+    if (totalnum == 0) {
+      break;
+    }
+    // If one of the cluster gets more than m_imageThreshold, divide
+    int divide = 0;
+    for (int i = 0; i < (int)m_timages.size(); ++i)
+      if (m_imageThreshold < (int)m_timages[i].size()) {
+        divide = 1;
+        break;
+      }
+    if (divide) {
+      break;
+    }
+    
+    setClusters();
+  }
+
+  for (int i = 0; i < (int)m_addnums.size(); ++i)
+    cout << m_addnums[i] << ' ';
+  cout << endl;
+  
+  int totalnum = 0;
+  for (int c = 0; c < (int)m_timages.size(); ++c)
+    totalnum += (int)m_timages[c].size();
+  int beforenum = 0;
+  for (int c = 0; c < m_cnum; ++c)
+    if (m_removed[c] == 0)
+      ++beforenum;
+  
+  cout << "Image nums: "
+       << m_cnum << " -> " << beforenum <<  " -> " << totalnum << endl;
+}
+
+int Cbundle::addImages(void) {
+  m_thread = 0;
+  m_jobs.clear();
+  
+  // we think about sfm points that belong to images that have not been satisfied
+  vector<char> flags;
+  flags.resize((int)m_coords.size(), 0);
+  for (int c = 0; c < m_cnum; ++c) {
+    if (m_lacks[c] <= 0)
+      continue;
+    else {
+      for (int i = 0; i < (int)m_vpoints[c].size(); ++i) {
+        const int pid = m_vpoints[c][i];
+        if (m_sfms2[pid].m_satisfied == 2)
+          flags[pid] = 1;
+      }
+    }
+  }
+
+#pragma omp parallel for
+  for (int p = 0; p < (int)m_sfms2.size(); ++p) {
+    if (flags[p] == 0)
+      continue;
+    m_sfms2[p].m_adds.clear();
+    
+    const int cluster = m_sfms2[p].m_cluster;
+    // Try to add an image to m_timages[cluster]
+    vector<int> cimages;
+    set_intersection(m_timages[cluster].begin(), m_timages[cluster].end(),
+                     m_visibles[p].begin(), m_visibles[p].end(),
+                     back_inserter(cimages));
+    sort(cimages.begin(), cimages.end());
+    
+    for (int i = 0; i < (int)m_visibles[p].size(); ++i) {
+      const int image = m_visibles[p][i];
+      
+      if (binary_search(cimages.begin(), cimages.end(), image))
+        continue;
+      
+      vector<int> vitmp = cimages;
+      vitmp.push_back(image);
+      const float newscore = computeScore2(m_coords[p], vitmp);
+      if (newscore <= m_sfms2[p].m_score)
+        continue;
+
+      m_sfms2[p].m_adds.push_back(Sadd(image,
+                                       (newscore - m_sfms2[p].m_score) /
+                                       m_sfms2[p].m_scoreThreshold));
+    }
+  }
+  
+  // Accumulate information from SfM points. For each cluster.
+  // For each cluster, for each image, sum of gains.
+  vector<map<int, float> > cands;
+  cands.resize((int)m_timages.size());
+
+  for (int p = 0; p < (int)m_sfms2.size(); ++p) {
+    if (flags[p] == 0)
+      continue;
+    
+    const int cluster = m_sfms2[p].m_cluster;
+    
+    for (int i = 0; i < (int)m_sfms2[p].m_adds.size(); ++i) {
+      Sadd& stmp = m_sfms2[p].m_adds[i];
+      cands[cluster][stmp.m_image] += stmp.m_gain;
+    }
+  }
+
+  return addImagesSub(cands);
+}
+
+int Cbundle::addImagesSub(const std::vector<std::map<int, float> >& cands) {
+  // Vec3f (gain, cluster, image). Start adding starting from the good
+  // one, block neighboring images.
+  vector<Vec3f> cands2;
+  for (int i = 0; i < (int)m_timages.size(); ++i) {
+    map<int, float>::const_iterator mbegin = cands[i].begin();
+    map<int, float>::const_iterator mend = cands[i].end();
+    
+    while (mbegin != mend) {
+      cands2.push_back(Vec3f(-mbegin->second, i, mbegin->first));
+      ++mbegin;
+    }
+  }
+
+  if (cands2.empty())
+    return 0;
+
+  sort(cands2.begin(), cands2.end(), Svec3cmp<float>());
+  
+  vector<char> blocked;
+  blocked.resize(m_cnum, 0);
+  vector<int> addnum;
+  addnum.resize((int)m_timages.size(), 0);
+
+  // A bit of tuning is possible here. For the paper, we used 0.7f,
+  //but 0.9f seems to produce better results in general.  const float
+  //gainThreshold = -cands2[0][0] * 0.7f;
+  const float gainThreshold = -cands2[0][0] * 0.9f;
+
+  // use rato threshold for blocked one. if not blocked, just keep on
+  // going.
+  for (int i = 0; i < (int)cands2.size(); ++i) {
+    const float gain = -cands2[i][0];
+    if (gain < gainThreshold)
+      break;
+    
+    const int cluster = (int)cands2[i][1];
+    const int image = (int)cands2[i][2];
+    
+    if (blocked[image])
+      continue;
+
+    // Add
+    ++addnum[cluster];
+    ++m_addnums[cluster];
+    blocked[image] = 1;
+    for (int i = 0; i < (int)m_neighbors[image].size(); ++i)
+      blocked[m_neighbors[image][i]] = 1;
+
+    m_timages[cluster].push_back(image);
+    // m_score, m_cluster, m_satisfied, m_lacks are updated in
+    // setClusters.  So, no need to update anything.
+  }
+  
+  for (int i = 0; i < (int)m_timages.size(); ++i)
+    sort(m_timages[i].begin(), m_timages[i].end());
+  
+  return accumulate(addnum.begin(), addnum.end(), 0);
+}
+
+int Cbundle::totalNum(void) const{
+  int totalnum = 0;
+  for (int c = 0; c < (int)m_timages.size(); ++c)
+    totalnum += (int)m_timages[c].size();
+  return totalnum;
+}
+
+int Cbundle::my_isIntersect(const std::vector<int>& lhs,
+                            const std::vector<int>& rhs) {
+  vector<int>::const_iterator b0 = lhs.begin();
+  vector<int>::const_iterator e0 = lhs.end();
+
+  vector<int>::const_iterator b1 = rhs.begin();
+  vector<int>::const_iterator e1 = rhs.end();
+
+  while (1) {
+    if (b0 == e0)
+      return 0;  
+    if (b1 == e1)
+      return 0;
+
+    if (*b0 == *b1) {
+      return 1;
+    }
+    else if (*b0 < *b1)
+      ++b0;
+    else
+      ++b1;
+  }
+}
+
+void Cbundle::mymerge(const std::vector<int>& lhs,
+                      const std::vector<int>& rhs,
+                      std::vector<int>& output) {
+  vector<int>::const_iterator b0 = lhs.begin();
+  vector<int>::const_iterator e0 = lhs.end();
+
+  vector<int>::const_iterator b1 = rhs.begin();
+  vector<int>::const_iterator e1 = rhs.end();
+
+  while (1) {
+    if (b0 == e0) {
+      output.insert(output.end(), b1, e1);
+      break;
+    }
+    if (b1 == e1) {
+      output.insert(output.end(), b0, e0);
+      break;
+    }
+
+    if (*b0 == *b1) {
+      output.push_back(*b0);
+      ++b0;
+      ++b1;
+    }
+    else if (*b0 < *b1) {
+      output.push_back(*b0);
+      ++b0;
+    }
+    else {
+      output.push_back(*b1);
+      ++b1;
+    }
+  }
+}
+
+// Calculates log2 of number.
+template <typename T>
+T Log2( T n )
+{
+    // log(n)/log(2) is log2.
+    return log( n ) / log( T(2) );
+}
+
+void Cbundle::setWidthsHeightsLevels(void) {
+  m_widths.resize(m_cnum);
+  m_heights.resize(m_cnum);
+  m_levels.resize(m_cnum);
+
+  // Determine level. SfM was done on 2M pixels.
+  for (int c = 0; c < m_cnum; ++c) {
+    m_levels[c] = m_dlevel;
+
+    m_widths[c] = m_pss.m_photos[c].getWidth(m_levels[c]);
+    m_heights[c] = m_pss.m_photos[c].getHeight(m_levels[c]);
+  }
+}
+
+ 
+float Cbundle::computeScore2(const Vec4f& coord,
+                             const std::vector<int>& images) const {
+  vector<int> uimages;
+  return computeScore2(coord, images, uimages);
+}
+
+float Cbundle::computeScore2(const Vec4f& coord,
+                             const std::vector<int>& images,
+                             std::vector<int>& uimages) const {
+  const int inum = (int)images.size();
+  uimages.clear();
+  if (inum < 2)
+    return -1.0f;
+  
+  vector<Vec4f> rays;    rays.resize(inum);
+  vector<float> scales;  scales.resize(inum);
+  
+  for (int r = 0; r < inum; ++r) {
+    rays[r] = m_pss.m_photos[images[r]].m_center - coord;
+    unitize(rays[r]);
+    
+    scales[r] = 1.0f / m_pss.m_photos[images[r]].getScale(coord, 0);
+  }
+
+  // Find the best pair of images
+  Vec2i bestpair;
+  float bestscore = -INT_MAX/2;
+  for (int i = 0; i < inum; ++i)
+    for (int j = i+1; j < inum; ++j) {
+      const float ftmp =
+        angleScore(rays[i], rays[j]) * scales[i] * scales[j];
+      if (bestscore < ftmp) {
+        bestscore = ftmp;
+        bestpair = Vec2i(i, j);
+      }
+    }
+
+  vector<int> inout;  inout.resize(inum, 1);
+  inout[bestpair[0]] = 0;         inout[bestpair[1]] = 0;
+  uimages.push_back(bestpair[0]);   uimages.push_back(bestpair[1]);
+
+  for (int t = 2; t < min(m_tau, inum); ++t) {
+    // Add the best new image
+    int ansid = -1;
+    float ansscore = -INT_MAX/2;
+    for (int j = 0; j < inum; ++j) {
+      if (inout[j] == 0)
+        continue;
+
+      float score = 0.0f;
+      for (int k = 0; k < (int)uimages.size(); ++k) {
+        const int iid = uimages[k];
+        score += angleScore(rays[j], rays[iid]) * scales[j] * scales[iid];
+      }
+      if (ansscore < score) {
+        ansscore = score;
+        ansid = j;
+      }
+    }
+    
+    if (ansid == -1) {
+      cerr << "Impossible 2 in compureScore" << endl;      exit (1);
+    }
+    
+    inout[ansid] = 0;
+    uimages.push_back(ansid);
+    bestscore += ansscore;
+  }
+
+  for (int i = 0; i < (int)uimages.size(); ++i)
+    uimages[i] = images[uimages[i]];
+  
+  return bestscore;
+}
+
+float Cbundle::computeScore2(const Vec4f& coord,
+                             const std::vector<int>& images,
+                             const int index) const {
+  vector<int> uimages;
+  const int inum = (int)images.size();
+  if (inum < 2)
+    return -1.0f;
+  
+  vector<Vec4f> rays;    rays.resize(inum);
+  vector<float> scales;  scales.resize(inum);
+  
+  for (int r = 0; r < inum; ++r) {
+    rays[r] = m_pss.m_photos[images[r]].m_center - coord;
+    unitize(rays[r]);
+    
+    scales[r] = 1.0f / m_pss.m_photos[images[r]].getScale(coord, 0);
+  }
+
+  float bestscore = 0.0f;
+  vector<int> inout;  inout.resize(inum, 1);
+  inout[index] = 0;
+  uimages.push_back(index);
+
+  for (int t = 1; t < min(m_tau, inum); ++t) {
+    // Add the best new image
+    int ansid = -1;
+    float ansscore = -INT_MAX/2;
+    for (int j = 0; j < inum; ++j) {
+      if (inout[j] == 0)
+        continue;
+
+      float score = 0.0f;
+      for (int k = 0; k < (int)uimages.size(); ++k) {
+        const int iid = uimages[k];
+        score += angleScore(rays[j], rays[iid]) * scales[j] * scales[iid];
+      }
+      if (ansscore < score) {
+        ansscore = score;
+        ansid = j;
+      }
+    }
+    
+    if (ansid == -1) {
+      cerr << "Impossible 2 in compureScore" << endl;      exit (1);
+    }
+    
+    inout[ansid] = 0;
+    uimages.push_back(ansid);
+    bestscore += ansscore;
+  }
+  return bestscore;
+}
+
+void Cbundle::writeVis(void) {
+  ofstream ofstr;
+  char buffer[1024];
+  sprintf(buffer, "%svis.dat", m_prefix.c_str());
+  
+  ofstr.open(buffer);
+  ofstr << "VISDATA" << endl;
+  ofstr << m_cnum << endl;
+
+  int numer = 0;
+  int denom = 0;
+  
+  for (int c = 0; c < m_cnum; ++c) {
+    if (m_removed[c])
+      ofstr << c << ' ' << 0 << endl;
+    else {
+      ofstr << c << ' ' << (int)m_neighbors[c].size() << "  ";
+      for (int i = 0; i < (int)m_neighbors[c].size(); ++i)
+        ofstr << m_neighbors[c][i] << ' ';
+      ofstr << endl;
+
+      numer += (int)m_neighbors[c].size();
+      ++denom;
+    }
+  }
+  ofstr.close();
+
+
+  cout << numer / (float)denom << " images in vis on the average" << endl;
+  
+}
+
+void Cbundle::writeCameraCenters(void) {
+  for (int i = 0; i < (int)m_timages.size(); ++i) {
+    char buffer[1024];
+    sprintf(buffer, "%scenters-%04d.ply", m_prefix.c_str(), i);
+    
+    ofstream ofstr;
+    ofstr.open(buffer);
+    ofstr << "ply" << endl
+          << "format ascii 1.0" << endl
+          << "element vertex " << (int)m_timages[i].size() << endl
+          << "property float x" << endl
+          << "property float y" << endl
+          << "property float z" << endl
+          << "end_header" << endl;
+    for (int j = 0; j < (int)m_timages[i].size(); ++j) {
+      const int c = m_timages[i][j];
+      ofstr << m_pss.m_photos[c].m_center[0] << ' '
+            << m_pss.m_photos[c].m_center[1] << ' '
+            << m_pss.m_photos[c].m_center[2] << endl;
+    }
+    ofstr.close();
+  }
+
+  {
+    char buffer[1024];
+    sprintf(buffer, "%scenters-all.ply", m_prefix.c_str());
+    ofstream ofstr;
+    ofstr.open(buffer);
+    
+    ofstr << "ply" << endl
+          << "format ascii 1.0" << endl
+          << "element vertex " << m_cnum << endl
+          << "property float x" << endl
+          << "property float y" << endl
+          << "property float z" << endl
+          << "property uchar diffuse_red" << endl
+          << "property uchar diffuse_green" << endl
+          << "property uchar diffuse_blue" << endl
+          << "end_header" << endl;
+    for (int c = 0; c < m_cnum; ++c) {
+      ofstr << m_pss.m_photos[c].m_center[0] << ' '
+            << m_pss.m_photos[c].m_center[1] << ' '
+            << m_pss.m_photos[c].m_center[2] << ' ';
+      
+      if (m_removed[c])
+        ofstr << "0 255 0" << endl;
+      else
+        ofstr << "255 0 255" << endl;
+    }
+    ofstr.close();
+  }
+}
+
+void Cbundle::writeGroups(void) {
+  char buffer[1024];
+  sprintf(buffer, "%sske.dat", m_prefix.c_str());
+  ofstream ofstr;
+  ofstr.open(buffer);
+  ofstr << "SKE" << endl
+        << m_cnum << ' ' << (int)m_timages.size() << endl;
+  
+  for (int c = 0; c < (int)m_timages.size(); ++c) {
+    ofstr << (int)m_timages[c].size() << ' ' << (int)m_oimages[c].size() << endl;
+
+    for (int i = 0; i < (int)m_timages[c].size(); ++i)
+      ofstr << m_timages[c][i] << ' ';
+    ofstr << endl;
+    for (int i = 0; i < (int)m_oimages[c].size(); ++i)
+      ofstr << m_oimages[c][i] << ' ';
+    ofstr << endl;
+  }
+}
+
+void Cbundle::startTimer(void) {
+  time(&m_tv); 
+}
+
+time_t Cbundle::curTimer(void) {
+  time(&m_curtime); 
+  return m_tv - m_curtime;
+}

From 6afab443e6aa7e4ae00fa5032086d98464ac2007 Mon Sep 17 00:00:00 2001
From: martin <cev_madde@msn.com>
Date: Thu, 18 Feb 2016 00:30:43 +0100
Subject: [PATCH 19/19] cerr breaks at this point

---
 program/base/cmvs/bundle.cc~ | 1476 ----------------------------------
 1 file changed, 1476 deletions(-)
 delete mode 100644 program/base/cmvs/bundle.cc~

diff --git a/program/base/cmvs/bundle.cc~ b/program/base/cmvs/bundle.cc~
deleted file mode 100644
index d0be6d08..00000000
--- a/program/base/cmvs/bundle.cc~
+++ /dev/null
@@ -1,1476 +0,0 @@
-#include <fstream>
-#include <iterator>
-#include <numeric> //PM
-
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-
-#include "graclus.h"
-#include "bundle.h"
-#define _USE_MATH_DEFINES
-#include <math.h>
-
-extern "C"
-{
-int boundary_points;
-int spectral_initialization = 0;
-int cutType;
-int memory_saving;
-};
-
-using namespace std;
-using namespace boost;
-using namespace CMVS;
-
-Cbundle::Cbundle(void) {
-  m_CPU = 8;
-  m_junit = 100;
-  mtx_init(&m_lock, mtx_plain | mtx_recursive);
-  m_debug = 0;
-  m_puf = NULL;
-  m_ptree = NULL;
-}
-
-Cbundle::~Cbundle() {
-}
-
-void Cbundle::prep(const std::string prefix, const int imageThreshold,
-                   const int tau, const float scoreRatioThreshold,
-                   const float coverageThreshold,
-                   const int pnumThreshold, const int CPU) {
-  if (pnumThreshold != 0) {
-    cerr << "Should use pnumThreshold = 0" << endl;
-    exit (1);
-  }
-  
-  m_prefix = prefix;
-  m_imageThreshold = imageThreshold;
-  m_tau = tau;
-  m_scoreRatioThreshold = scoreRatioThreshold;
-  m_coverageThreshold = coverageThreshold;
-  m_pnumThreshold = pnumThreshold;
-  m_CPU = CPU;
-  
-  m_linkThreshold = 2.0f;
-
-  m_dscale = 1 / 100.0f;
-  m_dscale2 = 1.0f;
-  
-  char buffer[1024];
-  sprintf(buffer, "%sbundle.rd.out", prefix.c_str());
-  cout << "Reading bundle..." << flush;
-  readBundle(buffer);
-  cout << endl;
-
-  vector<int> images;
-  for (int c = 0; c < m_cnum; ++c)
-    images.push_back(c);
-
-  m_dlevel = 7;
-  m_maxLevel = 12;
-  m_pss.init(images, prefix, m_maxLevel + 1, 5, 0);
-  
-  cout << "Set widths/heights..." << flush;
-  setWidthsHeightsLevels();
-  cout << "done" << flush;
-}
-
-void Cbundle::prep2(void) {
-  // Used in mergeSfMP now.
-  {
-    m_pweights.resize((int)m_coords.size(), 1);
-    m_sfms2.resize((int)m_coords.size());
-    startTimer();
-    setScoreThresholds();
-    cout << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
-    startTimer();
-    cout << "slimNeighborsSetLinks..." << flush;
-    slimNeighborsSetLinks();
-    cout << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
-  }
-  
-  // Improve visibility by using texture analysis
-  startTimer();  
-  cout << "mergeSFM..." << flush;
-  mergeSfMP();
-  cout << '\t' << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
-
-  m_sfms2.clear();
-  m_sfms2.resize((int)m_coords.size());
-  cout << "setScoreThresholds..." << flush;
-  startTimer();
-  setScoreThresholds();
-  cout << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
-
-  // Remove redundant images first
-  cout << "sRemoveImages... " << flush;
-  startTimer();
-
-  sRemoveImages();
-  cout << '\t' << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
-
-  // use m_removed to change m_visibles and update m_neighbors
-  startTimer();
-  resetVisibles();
-  setNeighbors();
-  cout << "slimNeighborsSetLinks..." << flush;
-  slimNeighborsSetLinks();
-  cout << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
-  
-  // Init m_timages by mutually exclusive clustering
-  setTimages();
-  m_oimages.resize((int)m_timages.size());  
-}
-
-void Cbundle::run(const std::string prefix, const int imageThreshold,
-                  const int tau, const float scoreRatioThreshold,
-                  const float coverageThreshold,
-                  const int pnumThreshold, const int CPU) {
-  startTimer();
-  
-  prep(prefix, imageThreshold, tau, scoreRatioThreshold,
-       coverageThreshold, pnumThreshold, CPU);
-
-  cout << '\t' << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
-
-  prep2();
-  
-  // Assumed variables that must be set properly here
-  cout << "Adding images: " << endl;
-  startTimer();
-  // Add images
-  // Repeat until all the clusters become at most m_imageThreshold.
-  while (1) {
-    addImagesP();
-    
-    int change = 0;
-    vector<vector<int> > newtimages;
-    cout << "Divide: " << flush;
-    for (int i = 0; i < (int)m_timages.size(); ++i) {
-      if ((int)m_timages[i].size() <= m_imageThreshold) {        
-        newtimages.push_back(m_timages[i]);
-        continue;
-      }
-      else {
-        cout << i << ' ';
-        
-        change = 1;
-        // divide
-        vector<vector<int> > vvi;
-        divideImages(m_timages[i], vvi);
-        for (int j = 0; j < (int)vvi.size(); ++j)
-          newtimages.push_back(vvi[j]);
-      }
-    }
-
-    cout << endl;
-    
-    m_timages.swap(newtimages);
-    if (change == 0)
-      break;
-  }
-  cout << "done\t" << curTimer()/CLOCKS_PER_SEC << " secs" << endl;
-  
-  m_oimages.resize((int)m_timages.size());
-
-  writeCameraCenters();
-  
-  // Output results
-  writeVis();
-  writeGroups();
-}
-
-float Cbundle::computeLink(const int image0, const int image1) {
-  vector<int> common;
-  set_intersection(m_vpoints[image0].begin(), m_vpoints[image0].end(),
-                   m_vpoints[image1].begin(), m_vpoints[image1].end(),
-                   back_inserter(common));
-
-  float score = 0.0f;
-  for (int i = 0; i < (int)common.size(); ++i) {
-    const int pid = common[i];
-    vector<int> vtmp;
-    vtmp.push_back(image0);    vtmp.push_back(image1);
-    const float ftmp = computeScore2(m_coords[pid], vtmp);
-    if (m_sfms2[pid].m_score != 0.0f)
-      score += m_pweights[pid] *
-        ftmp / (m_sfms2[pid].m_scoreThreshold / m_scoreRatioThreshold);
-  }
-
-  return score;
-}
-
-void Cbundle::slimNeighborsSetLinks(void) {
-  const int maxneighbor = 30;
-  m_links.clear();
-  m_links.resize(m_cnum);
-
-#pragma omp parallel for
-  for (int image = 0; image < m_cnum; ++image) {
-    m_links[image].resize((int)m_neighbors[image].size());
-    
-    for (int i = 0; i < (int)m_neighbors[image].size(); ++i)
-      m_links[image][i] = computeLink(image, m_neighbors[image][i]);
-
-    if ((int)m_neighbors[image].size() < 2)
-      continue;
-
-    vector<int> newneighbors;
-    vector<float> newlinks;
-
-    vector<Vec2f> vv;
-    for (int i = 0; i < (int)m_neighbors[image].size(); ++i)
-      vv.push_back(Vec2(-m_links[image][i], m_neighbors[image][i]));
-    sort(vv.begin(), vv.end(), Svec2cmp<float>());
-
-    const int itmp = min(maxneighbor, (int)vv.size());
-    for (int i = 0; i < itmp; ++i) {
-      newneighbors.push_back((int)vv[i][1]);
-      newlinks.push_back(-vv[i][0]);
-    }
-    
-    m_neighbors[image].swap(newneighbors);
-    m_links[image].swap(newlinks);
-  }
-}
-
-void Cbundle::setScoreThresholds(void) {
-#pragma omp parallel for
-  for (int p = 0; p < (int)m_coords.size(); ++p) {
-    m_sfms2[p].m_scoreThreshold =
-      computeScore2(m_coords[p], m_visibles[p], m_sfms2[p].m_uimages)
-      * m_scoreRatioThreshold;
-  }
-}
-
-void Cbundle::sRemoveImages(void) {
-  m_removed.clear();
-  m_removed.resize(m_cnum, 0);
-    
-  m_allows.resize(m_cnum);
-  for (int c = 0; c < m_cnum; ++c)
-    m_allows[c] = (int)ceil((int)m_vpoints[c].size() *
-                            (1.0f - m_coverageThreshold));
-  
-  // Sort all the images in an increasing order of resolution
-  vector<Vec2i> vvi;
-  for (int c = 0; c < m_cnum; ++c) {
-    const int res =
-      m_pss.m_photos[c].getWidth(0) * m_pss.m_photos[c].getHeight(0);
-    vvi.push_back(Vec2i(res, c));
-  }
-  sort(vvi.begin(), vvi.end(), Svec2cmp<int>());
-  
-  const int tenth = max(1, m_cnum / 10);
-  for (int i = 0; i < (int)vvi.size(); ++i) {
-    if (i % tenth == 0)
-      cout << '*' << flush;
-    const int image = vvi[i][1];
-    checkImage(image);
-  }
-  cout << endl;
-    
-  cout << "Kept: ";
-  int kept = 0;
-  for (int c = 0; c < m_cnum; ++c)
-    if (m_removed[c] == 0) {
-      ++kept;
-      cout << c << ' ';
-    }
-  cout << endl << endl;
-  
-  cout << "Removed: ";
-  for (int c = 0; c < m_cnum; ++c)
-    if (m_removed[c]) {
-      cout << c << ' ';
-    }
-  cout << endl;
-  cout << "sRemoveImages: " << m_cnum << " -> " << kept << flush;
-}
-
-void Cbundle::resetVisibles(void) {
-  // reset m_visibles. remove "removed images" from the list.
-  for (int p = 0; p < (int)m_visibles.size(); ++p) {
-    vector<int> newimages;
-
-    setNewImages(p, -1, newimages);
-    m_visibles[p].swap(newimages);
-  }
-}
-
-void Cbundle::setNewImages(const int pid, const int rimage,
-                           std::vector<int>& newimages) {
-  for (int i = 0; i < (int)m_visibles[pid].size(); ++i) {
-    const int image = m_visibles[pid][i];
-    if (m_removed[image] || image == rimage)
-      continue;
-    newimages.push_back(image);
-  }
-}
-
-void Cbundle::checkImage(const int image) {
-  // For each SfM, check if it will be unsatisfied by removing image.
-  //  0: unsatisfy
-  //  1: satisfy->satisfy
-  //  2: satisfy->unsatisfy
-  m_statsT.resize((int)m_vpoints[image].size());
-  m_jobs.clear();
-  
-#pragma omp parallel for
-  for (int p = 0; p < (int)m_statsT.size(); ++p) {
-    const int pid = m_vpoints[image][p];
-    if (m_sfms2[pid].m_satisfied == 0)
-      m_statsT[p] = 0;
-    else {
-      m_statsT[p] = 1;
-
-      // if m_uimages of p is not removed, if image is not in
-      // m_uimages, the point is still satisfied.
-      int valid = 1;
-      int inside = 0;
-
-      for (int i = 0; i < (int)m_sfms2[pid].m_uimages.size(); ++i) {
-        const int itmp = m_sfms2[pid].m_uimages[i];
-        if (itmp == image)
-          inside = 1;
-        if (m_removed[itmp]) {
-          valid = 0;
-          break;
-        }
-      }
-      
-      if (valid == 1 && inside == 0)
-        continue;
-      
-      vector<int> newimages;
-      setNewImages(pid, image, newimages);
-      const float cscore = computeScore2(m_coords[pid], newimages);
-      if (cscore < m_sfms2[pid].m_scoreThreshold)
-        m_statsT[p] = 2;
-    }
-  }
-  
-  // For each image, how many SFM points are removed if you remove
-  // "image"
-  vector<int> decrements;
-  decrements.resize(m_cnum, 0);
-  
-  // Look at points with m_stastT[p] = 2, to see if m_allows are still
-  // above thresholds.
-  for (int p = 0; p < (int)m_statsT.size(); ++p) {
-    if (m_statsT[p] == 2) {
-      const int pid = m_vpoints[image][p];
-      for (int i = 0; i < (int)m_visibles[pid].size(); ++i) {
-        const int itmp = m_visibles[pid][i];
-        ++decrements[itmp];
-      }
-    }
-  }
-
-  // If m_allows can cover decrements, go ahead.
-  int rflag = 1;
-  for (int c = 0; c < m_cnum; ++c)
-    if (m_allows[c] < decrements[c]) {
-      rflag = 0;
-      break;
-    }
-
-  // remove
-  if (rflag) {
-    m_removed[image] = 1;
-    for (int p = 0; p < (int)m_statsT.size(); ++p)
-      if (m_statsT[p] == 2)
-        m_sfms2[m_vpoints[image][p]].m_satisfied = 0;
-
-    for (int c = 0; c < m_cnum; ++c)
-      m_allows[c] -= decrements[c];
-
-    // Update uimages for points that are still satisfied and still
-    // contains image in m_uimages
-#pragma omp parallel for
-    for (int p = 0; p < (int)m_statsT.size(); ++p) {
-      const int pid = m_vpoints[image][p];
-      if (m_statsT[p] == 1) {
-        int contain = 0;
-        for (int i = 0; i < (int)m_sfms2[pid].m_uimages.size(); ++i)
-          if (m_sfms2[pid].m_uimages[i] == image) {
-            contain = 1;
-            break;
-          }
-
-        if (contain) {
-          vector<int> newimages;
-          setNewImages(pid, -1, newimages);
-          const float cscore =
-            computeScore2(m_coords[pid], newimages, m_sfms2[pid].m_uimages);
-          if (cscore < m_sfms2[pid].m_scoreThreshold)
-            m_sfms2[pid].m_satisfied = 0;
-        }
-      }
-    }
-  }
-}
-
-void Cbundle::setNeighbors(void) {
-  m_neighbors.clear();
-  m_neighbors.resize(m_cnum);
-#pragma omp parallel for
-  for (int image = 0; image < m_cnum; ++image) {
-    vector<int> narray;
-    narray.resize(m_cnum, 0);
-    
-    for (int p = 0; p < (int)m_visibles.size(); ++p) {
-      if (!binary_search(m_visibles[p].begin(), m_visibles[p].end(),
-                         image))
-        continue;
-      
-      for (int i = 0; i < (int)m_visibles[p].size(); ++i)
-        narray[m_visibles[p][i]] = 1;
-    }
-    
-    for (int i = 0; i < m_cnum; ++i)
-      if (narray[i] && i != image)
-        m_neighbors[image].push_back(i);
-  } 
-}
-
-void Cbundle::setTimages(void) {
-  vector<int> lhs;
-  for (int c = 0; c < m_cnum; ++c)
-    if (m_removed[c] == 0)
-      lhs.push_back(c);
-
-  m_timages.clear();
-
-  if ((int)lhs.size() <= m_imageThreshold)
-    m_timages.push_back(lhs);
-  else
-    divideImages(lhs, m_timages);
-  
-  cout << endl << "Cluster sizes: " << endl;
-  for (int i = 0; i < (int)m_timages.size(); ++i)
-    cout << (int)m_timages[i].size() << ' ';
-  cout << endl;
-}
-
-void Cbundle::divideImages(const std::vector<int>& lhs,
-                           std::vector<std::vector<int> >& rhs) {
-  const float iratio = 125.0f / 150.0f;
-  // Candidates
-  list<vector<int> > candidates;
-  // initialize first cluster
-  candidates.push_back(vector<int>());
-  vector<int>& vitmp = candidates.back();
-  for (int i = 0; i < (int)lhs.size(); ++i)
-    vitmp.push_back(lhs[i]);
-  
-  // Iterate
-  while (1) {
-    if (candidates.empty())
-      break;
-    vector<int> cand = candidates.front();
-    candidates.pop_front();
-    if ((int)cand.size() <= m_imageThreshold * iratio) {
-      rhs.push_back(cand);
-      continue;
-    }
-    // Divide into 2 groups
-    vector<idxtype> xadj, adjncy, adjwgt, part;
-    const int nparts = 2;
-    const int cutType = 0;
-    // Define graphs
-    for (int i = 0; i < (int)cand.size(); ++i) {
-      xadj.push_back((int)adjncy.size());
-
-      for (int j = 0; j < (int)cand.size(); ++j) {
-        if (i == j)
-          continue;
-        // Check if cand[i] and cand[j] are connected
-        const int cid0 = cand[i];
-        const int cid1 = cand[j];
-        vector<int>::const_iterator nite = 
-          find(m_neighbors[cid0].begin(), m_neighbors[cid0].end(), cid1);
-        
-        if (nite != m_neighbors[cid0].end()) {
-          adjncy.push_back(j);
-          const int offset = nite - m_neighbors[cid0].begin();
-          const int weight =
-            min(5000, (int)floor(10.0f * m_links[cid0][offset] + 0.5f));
-          adjwgt.push_back(weight);
-        }
-      }
-    }
-    xadj.push_back((int)adjncy.size());
-
-    Cgraclus::runE(xadj, adjncy, adjwgt, nparts, cutType, part);
-    
-    // Divide into 2 groups
-    vector<int> cand1, cand2;
-    for (int i = 0; i < (int)part.size(); ++i) {
-      if (part[i] == 0)
-        cand1.push_back(cand[i]);
-      else
-        cand2.push_back(cand[i]);
-    }
-    
-    if (cand1.empty() || cand2.empty()) {
-      cerr << "Error. Normalized cuts produced an empty cluster: "
-           << (int)part.size() << " -> "
-           << (int)cand1.size() << ' '
-           << (int)cand2.size() << endl;
-      exit (1);
-    }
-    
-    if ((int)cand1.size() <= m_imageThreshold * iratio)
-      rhs.push_back(cand1);
-    else {
-      candidates.push_back(vector<int>());
-      (candidates.back()).swap(cand1);
-    }
-    if ((int)cand2.size() <= m_imageThreshold * iratio)
-      rhs.push_back(cand2);
-    else {
-      candidates.push_back(vector<int>());
-      (candidates.back()).swap(cand2);
-    }
-  }
-}
-
-void Cbundle::readBundle(const std::string file) {
-  // For each valid image, a list of connected images, and the second
-  // value is the number of common points.
-  ifstream ifstr;  ifstr.open(file.c_str());
-  if (!ifstr.is_open()) {
-    cerr << "Bundle file not found: " << file << endl;
-    exit (1);
-  }
-  while (1) {
-    unsigned char uctmp;
-    ifstr.read((char*)&uctmp, sizeof(unsigned char));
-    ifstr.putback(uctmp);
-    if (uctmp == '#') {
-      char buffer[1024];      ifstr.getline(buffer, 1024);
-    }
-    else
-      break;
-  }
-  int cnum, pnum;
-  ifstr >> cnum >> pnum;
-  vector<int> ids;        ids.resize(cnum);
-  cout << cnum << " cameras -- " << pnum << " points in bundle file" << endl;
-  m_cnum = 0;
-  for (int c = 0; c < cnum; ++c) {
-    ids[c] = -1;
-    float params[15];
-    for (int i = 0; i < 15; ++i)
-      ifstr >> params[i];
-    if (params[0] != 0.0f)
-      ids[c] = m_cnum++;
-  }
-
-  m_vpoints.resize(m_cnum);
-
-  m_coords.clear();  m_visibles.clear();
-  m_colors.clear();
-  m_pnum = 0;
-  const int tenth = max(1, pnum / 10);
-  for (int p = 0; p < pnum; ++p) {
-    if (p % tenth == 0)
-      cout << '*' << flush;
-    int num;    Vec3f color;
-    Vec4f coord;
-    ifstr >> coord[0] >> coord[1] >> coord[2]
-          >> color >> num;
-    coord[3] = 1.0f;
-    
-    vector<int> visibles;
-    for (int i = 0; i < num; ++i) {
-      int itmp;      ifstr >> itmp;
-      if (cnum <= itmp) {
-        double dtmp;
-        ifstr >> itmp >> dtmp >> dtmp;
-        continue;
-      }
-      if (ids[itmp] == -1) {
-        cerr << "impossible " << itmp << ' ' << ids[itmp] << endl;
-        exit (1);
-      }
-      visibles.push_back(ids[itmp]);
-
-      // Based on the bundler version, the number of parameters here
-      // are either 1 or 3. Currently, it seems to be 3.
-      double dtmp;
-      ifstr >> itmp >> dtmp >> dtmp;
-      //ifstr >> dtmp;
-    }
-    if ((int)visibles.size() < 2)
-      continue;
-
-    sort(visibles.begin(), visibles.end());
-    
-    for (int i = 0; i < (int)visibles.size(); ++i) {
-       m_vpoints[visibles[i]].push_back(m_pnum);
-    }
-    
-    m_visibles.push_back(visibles);
-    m_coords.push_back(coord);
-    m_colors.push_back(color);
-    ++m_pnum;
-  }
-  ifstr.close();
-  setNeighbors();
-  
-  cout << endl << m_cnum << " cameras -- " << m_pnum << " points" << flush;
-}
-
-void Cbundle::findPNeighbors(sfcnn<const float*, 3, float>& tree,
-                             const int pid, std::vector<int>& pneighbors) {
-  const float thresh = m_dscale2 * m_dscale2 *
-    m_minScales[pid] * m_minScales[pid];
-
-  vector<long unsigned int> ids;
-  vector<double> dists;
-  int k = min((int)m_coords.size(), 400);
-
-  while (1) {
-    ids.clear();
-    dists.clear();
-    tree.ksearch(&m_coords[pid][0], k, ids, dists);
-
-    if (thresh < dists[k - 1])
-      break;
-    if (k == (int)m_coords.size())
-      break;
-    
-    k = min((int)m_coords.size(), k * 2);
-  }
-
-  for (int i = 0; i < k; ++i) {
-    if (thresh < dists[i])
-      break;
-
-    // If distance from pid2 is also within threshold ok
-    const float thresh2 = m_dscale2 * m_dscale2 *
-      m_minScales[ids[i]] * m_minScales[ids[i]];
-    if (thresh2 < dists[i])
-      continue;
-    
-    if (ids[i] != (unsigned int)pid)
-      pneighbors.push_back(ids[i]);
-  }
-}
-
-void Cbundle::mergeSfMPThread(void) {
-  const int tenth = (int)m_coords.size() / 10;
-  while (1) {
-    int pid = -1;
-    mtx_lock(&m_lock);
-    if (!m_jobs.empty()) {
-      pid = m_jobs.front();
-      m_jobs.pop_front();
-    }
-    if (pid != -1 && m_merged[pid])
-      pid = -2;
-    if (m_count % tenth == 0)
-      cout << '*' << flush;
-    ++m_count;
-    mtx_unlock(&m_lock);
-    if (pid == -2)
-      continue;
-    if (pid == -1)
-      break;
-
-    // Process, and check later if merged flag is on
-    vector<int> pneighbors;
-    findPNeighbors(*m_ptree, pid, pneighbors);
-    const int psize = (int)pneighbors.size();
-    // visible images and its neighbors
-    vector<int> visibles = m_visibles[pid];
-    vector<int> vitmp;
-    for (int i = 0; i < (int)m_visibles[pid].size(); ++i) {
-      const int imagetmp = m_visibles[pid][i];
-      vitmp.clear();
-      mymerge(visibles, m_neighbors[imagetmp], vitmp);
-      vitmp.swap(visibles);
-    }
-    
-    vector<char> visflag(psize, 0);
-    // test visibility
-    for (int i = 0; i < psize; ++i) {
-      const int& pid2 = pneighbors[i];
-      if (my_isIntersect(visibles, m_visibles[pid2]))
-        visflag[i] = 1;
-    }
-    
-    // Now lock and try to register
-    mtx_lock(&m_lock);
-    // If the main one is removed, over... waste.
-    if (m_merged[pid] == 0) {
-      m_merged[pid] = 1;
-      for (int i = 0; i < psize; ++i) {
-        const int pid2 = pneighbors[i];
-        if (visflag[i] && m_merged[pid2] == 0) {
-          m_merged[pid2] = 1;
-          m_puf->union_set(pid, pid2);
-        }
-      }
-    }
-    mtx_unlock(&m_lock);
-  }
-}
-
-int Cbundle::mergeSfMPThreadTmp(void* arg) {
-  ((Cbundle*)arg)->mergeSfMPThread();
-  return 0;
-}
-
-void Cbundle::mergeSfMP(void) {
-  // Repeat certain times until no change
-  const int cpnum = (int)m_coords.size();
-  m_minScales.resize(cpnum);
-  for (int p = 0; p < cpnum; ++p) {
-    m_minScales[p] = INT_MAX/2;
-    for (int i = 0; i < (int)m_visibles[p].size(); ++i) {
-      const int image = m_visibles[p][i];
-      const float stmp =
-        m_pss.m_photos[image].getScale(m_coords[p], m_levels[image]);
-      m_minScales[p] = min(m_minScales[p], stmp);
-    }
-  }
-  
-  m_puf = new disjoint_sets_with_storage<>(cpnum);
-  for (int p = 0; p < cpnum; ++p)
-    m_puf->make_set(p);
-
-  {
-    // kdtree
-    vector<const float*> ppoints;
-    ppoints.resize((int)m_coords.size());
-    for (int i = 0; i < (int)m_coords.size(); ++i)
-      ppoints[i] = &(m_coords[i][0]);
-    m_ptree =
-      new sfcnn<const float*, 3, float>(&ppoints[0], (int)ppoints.size());
-    
-    m_merged.resize((int)m_coords.size(), 0);
-    m_jobs.clear();
-    vector<int> order;
-    order.resize(cpnum);
-    for (int p = 0; p < cpnum; ++p)
-      order[p] = p;
-    random_shuffle(order.begin(), order.end());
-    
-    for (int p = 0; p < cpnum; ++p)
-      m_jobs.push_back(order[p]);
-    
-    m_count = 0;
-    vector<thrd_t> threads(m_CPU);
-    for (int c = 0; c < m_CPU; ++c)
-      thrd_create(&threads[c], &mergeSfMPThreadTmp, (void*)this);
-    for (int c = 0; c < m_CPU; ++c)
-      thrd_join(threads[c], NULL);
-
-    int newpnum = 0;
-    // Mapping from m_pnum to new id for reps
-    vector<int> dict;
-    vector<int> reps;
-    dict.resize((int)m_coords.size(), -1);
-    for (int p = 0; p < (int)m_coords.size(); ++p) {
-      const int ptmp = m_puf->find_set(p);
-      if (p == ptmp) {
-        dict[p] = newpnum;
-        reps.push_back(p);
-        ++newpnum;
-      }
-    }
-  }
-  
-  // Based on m_puf, reset m_coords, m_coords, m_visibles, m_vpoints
-  cout << "resetPoints..." << flush;
-  resetPoints();
-  cout << "done" << endl;
-
-  delete m_puf;
-  m_puf = NULL;
-  delete m_ptree;
-  m_ptree = NULL;
-   
-  const int npnum = (int)m_coords.size();
-  cout << "Rep counts: " << cpnum << " -> " << npnum << "  " << flush;
-}
-
-// Based on m_puf, reset m_coords, m_coords, m_visibles, m_vpoints
-void Cbundle::resetPoints(void) {
-  vector<int> counts;
-  vector<int> smallestids;
-  counts.resize((int)m_coords.size(), 0);
-  smallestids.resize((int)m_coords.size(), INT_MAX/2);
-  for (int p = 0; p < (int)m_coords.size(); ++p) {
-    const int ptmp = m_puf->find_set(p);
-    smallestids[ptmp] = min(smallestids[ptmp], p);
-    ++counts[ptmp];
-  }
-  const int mthreshold = 2;
-  vector<Vec2i> vv2;
-  for (int p = 0; p < (int)m_coords.size(); ++p)
-    if (mthreshold <= counts[p])
-      vv2.push_back(Vec2i(smallestids[p], p));
-  sort(vv2.begin(), vv2.end(), Svec2cmp<int>());
-  
-  vector<int> newpids;
-  newpids.resize((int)m_coords.size(), -1);
-  int newpnum = (int)vv2.size();
-  for (int i = 0; i < newpnum; ++i)
-    newpids[vv2[i][1]] = i;
-  
-  vector<Vec4f> newcoords;  newcoords.resize(newpnum);
-  vector<vector<int> > newvisibles;  newvisibles.resize(newpnum);
-
-  for (int p = 0; p < (int)m_coords.size(); ++p) {
-    const int ptmp = m_puf->find_set(p);
-    if (counts[ptmp] < mthreshold)
-      continue;
-    
-    const int newpid = newpids[ptmp];
-    newcoords[newpid] += m_coords[p];
-    vector<int> vitmp;
-    mymerge(newvisibles[newpid], m_visibles[p], vitmp);
-    vitmp.swap(newvisibles[newpid]);
-  }  
-  
-  for (int i = 0; i < newpnum; ++i)
-    newcoords[i] = newcoords[i] / newcoords[i][3];
-  
-  // Update m_vpoints
-  m_coords.swap(newcoords);
-  m_visibles.swap(newvisibles);
-  
-  for (int c = 0; c < m_cnum; ++c)
-    m_vpoints[c].clear();  
-  
-  for (int p = 0; p < (int)m_coords.size(); ++p) {
-    for (int i = 0; i < (int)m_visibles[p].size(); ++i) {
-      const int itmp = m_visibles[p][i];
-      m_vpoints[itmp].push_back(p);
-    }
-  }
-
-  m_pweights.clear();
-  for (int i = 0; i < newpnum; ++i)
-    m_pweights.push_back(counts[vv2[i][1]]);
-}
-
-void Cbundle::setScoresClusters(void) {
-#pragma omp parallel for
-  for (int p = 0; p < (int)m_coords.size(); ++p) {
-    // if m_satisfied is 0, no hope (even all the images are in the
-    // cluster, not satisfied
-    if (m_sfms2[p].m_satisfied == 0)
-      continue;
-    
-    m_sfms2[p].m_satisfied = 2;
-    setCluster(p);
-  }
-}
-
-void Cbundle::setCluster(const int p) {
-  m_sfms2[p].m_score = -1.0f;
-  m_sfms2[p].m_cluster = -1;
-  for (int c = 0; c < (int)m_timages.size(); ++c) {
-    // Select images in cluster
-    vector<int> vitmp;
-    set_intersection(m_timages[c].begin(), m_timages[c].end(),
-                     m_visibles[p].begin(), m_visibles[p].end(),
-                     back_inserter(vitmp));
-    
-    const float stmp = computeScore2(m_coords[p], vitmp);
-    if (m_sfms2[p].m_score < stmp) {
-      m_sfms2[p].m_cluster = c;
-      m_sfms2[p].m_score = stmp;
-    }
-  }
-
-  // If no cluster contains 2 images
-  if (m_sfms2[p].m_cluster == -1) {
-    int find = 0;
-    for (int j = 0; j < (int)m_visibles[p].size(); ++j) {
-      if (find)
-        break;
-      for (int c = 0; c < (int)m_timages.size(); ++c)
-        if (binary_search(m_timages[c].begin(), m_timages[c].end(),
-                          m_visibles[p][j])) {
-          m_sfms2[p].m_cluster = c;
-          m_sfms2[p].m_score = 0.0f;
-          find = 1;
-          break;
-        }
-    }
-    // If none of the visibles images are
-    if (find == 0) {
-      cerr << "Impossible in setscoresclustersthread" << endl
-           << (int)m_visibles[p].size() << endl
-           << (int)m_sfms2[p].m_satisfied << endl;
-      exit (1);
-    }
-  }
-
-  if (m_sfms2[p].m_scoreThreshold <= m_sfms2[p].m_score) {
-    // SATISFIED
-    m_sfms2[p].m_satisfied = 1;
-
-    //  update m_lacks
-    mtx_lock(&m_lock);    
-    for (int i = 0; i < (int)m_visibles[p].size(); ++i) {
-      const int image = m_visibles[p][i];
-      --m_lacks[image];
-    }
-    mtx_unlock(&m_lock);
-  }
-}
-
-float Cbundle::angleScore(const Vec4f& ray0, const Vec4f& ray1) {
-  const static float lsigma = 5.0f * M_PI / 180.f;
-  const static float rsigma = 15.0f * M_PI / 180.0f;
-  const static float lsigma2 = 2.0f * lsigma * lsigma;
-  const static float rsigma2 = 2.0f * rsigma * rsigma;
-  const static float pivot = 20.0f * M_PI / 180.0f;
-
-  const float angle = acos(min(1.0f, ray0 * ray1));
-  const float diff = angle - pivot;
-
-  if (angle < pivot)
-    return exp(- diff * diff / lsigma2);
-  else
-    return exp(- diff * diff / rsigma2);
-}
-
-void Cbundle::setClusters(void) {
-#pragma omp parallel for
-  for (int p = 0; p < (int)m_sfms2.size(); ++p) {
-    if (m_sfms2[p].m_satisfied != 2)
-      continue;
-
-    setCluster(p);
-  }
-}
-
-void Cbundle::addImagesP(void) {
-  // set m_lacks (how many more sfm points need to be satisfied)
-  m_lacks.resize(m_cnum);
-  for (int c = 0; c < m_cnum; ++c) {
-    if (m_removed[c])
-      m_lacks[c] = 0;
-    else
-      m_lacks[c] =
-        (int)floor((int)m_vpoints[c].size() * m_coverageThreshold);
-  }
-
-  // Compute best score given all the images. Upper bound on the score
-  setScoresClusters();
-  
-  // Add images to clusters to make m_lacks at most 0 for all the
-  // images. In practice, for each cluster, identify corresponding sfm
-  // points that have not been satisfied. These points compute gains.
-  
-  // set m_adds for each sfm point, and add images
-
-  m_addnums.clear();
-  m_addnums.resize((int)m_timages.size(), 0);
-  
-  while (1) {
-    const int totalnum = addImages();
-    // Very end
-    if (totalnum == 0) {
-      break;
-    }
-    // If one of the cluster gets more than m_imageThreshold, divide
-    int divide = 0;
-    for (int i = 0; i < (int)m_timages.size(); ++i)
-      if (m_imageThreshold < (int)m_timages[i].size()) {
-        divide = 1;
-        break;
-      }
-    if (divide) {
-      break;
-    }
-    
-    setClusters();
-  }
-
-  for (int i = 0; i < (int)m_addnums.size(); ++i)
-    cout << m_addnums[i] << ' ';
-  cout << endl;
-  
-  int totalnum = 0;
-  for (int c = 0; c < (int)m_timages.size(); ++c)
-    totalnum += (int)m_timages[c].size();
-  int beforenum = 0;
-  for (int c = 0; c < m_cnum; ++c)
-    if (m_removed[c] == 0)
-      ++beforenum;
-  
-  cout << "Image nums: "
-       << m_cnum << " -> " << beforenum <<  " -> " << totalnum << endl;
-}
-
-int Cbundle::addImages(void) {
-  m_thread = 0;
-  m_jobs.clear();
-  
-  // we think about sfm points that belong to images that have not been satisfied
-  vector<char> flags;
-  flags.resize((int)m_coords.size(), 0);
-  for (int c = 0; c < m_cnum; ++c) {
-    if (m_lacks[c] <= 0)
-      continue;
-    else {
-      for (int i = 0; i < (int)m_vpoints[c].size(); ++i) {
-        const int pid = m_vpoints[c][i];
-        if (m_sfms2[pid].m_satisfied == 2)
-          flags[pid] = 1;
-      }
-    }
-  }
-
-#pragma omp parallel for
-  for (int p = 0; p < (int)m_sfms2.size(); ++p) {
-    if (flags[p] == 0)
-      continue;
-    m_sfms2[p].m_adds.clear();
-    
-    const int cluster = m_sfms2[p].m_cluster;
-    // Try to add an image to m_timages[cluster]
-    vector<int> cimages;
-    set_intersection(m_timages[cluster].begin(), m_timages[cluster].end(),
-                     m_visibles[p].begin(), m_visibles[p].end(),
-                     back_inserter(cimages));
-    sort(cimages.begin(), cimages.end());
-    
-    for (int i = 0; i < (int)m_visibles[p].size(); ++i) {
-      const int image = m_visibles[p][i];
-      
-      if (binary_search(cimages.begin(), cimages.end(), image))
-        continue;
-      
-      vector<int> vitmp = cimages;
-      vitmp.push_back(image);
-      const float newscore = computeScore2(m_coords[p], vitmp);
-      if (newscore <= m_sfms2[p].m_score)
-        continue;
-
-      m_sfms2[p].m_adds.push_back(Sadd(image,
-                                       (newscore - m_sfms2[p].m_score) /
-                                       m_sfms2[p].m_scoreThreshold));
-    }
-  }
-  
-  // Accumulate information from SfM points. For each cluster.
-  // For each cluster, for each image, sum of gains.
-  vector<map<int, float> > cands;
-  cands.resize((int)m_timages.size());
-
-  for (int p = 0; p < (int)m_sfms2.size(); ++p) {
-    if (flags[p] == 0)
-      continue;
-    
-    const int cluster = m_sfms2[p].m_cluster;
-    
-    for (int i = 0; i < (int)m_sfms2[p].m_adds.size(); ++i) {
-      Sadd& stmp = m_sfms2[p].m_adds[i];
-      cands[cluster][stmp.m_image] += stmp.m_gain;
-    }
-  }
-
-  return addImagesSub(cands);
-}
-
-int Cbundle::addImagesSub(const std::vector<std::map<int, float> >& cands) {
-  // Vec3f (gain, cluster, image). Start adding starting from the good
-  // one, block neighboring images.
-  vector<Vec3f> cands2;
-  for (int i = 0; i < (int)m_timages.size(); ++i) {
-    map<int, float>::const_iterator mbegin = cands[i].begin();
-    map<int, float>::const_iterator mend = cands[i].end();
-    
-    while (mbegin != mend) {
-      cands2.push_back(Vec3f(-mbegin->second, i, mbegin->first));
-      ++mbegin;
-    }
-  }
-
-  if (cands2.empty())
-    return 0;
-
-  sort(cands2.begin(), cands2.end(), Svec3cmp<float>());
-  
-  vector<char> blocked;
-  blocked.resize(m_cnum, 0);
-  vector<int> addnum;
-  addnum.resize((int)m_timages.size(), 0);
-
-  // A bit of tuning is possible here. For the paper, we used 0.7f,
-  //but 0.9f seems to produce better results in general.  const float
-  //gainThreshold = -cands2[0][0] * 0.7f;
-  const float gainThreshold = -cands2[0][0] * 0.9f;
-
-  // use rato threshold for blocked one. if not blocked, just keep on
-  // going.
-  for (int i = 0; i < (int)cands2.size(); ++i) {
-    const float gain = -cands2[i][0];
-    if (gain < gainThreshold)
-      break;
-    
-    const int cluster = (int)cands2[i][1];
-    const int image = (int)cands2[i][2];
-    
-    if (blocked[image])
-      continue;
-
-    // Add
-    ++addnum[cluster];
-    ++m_addnums[cluster];
-    blocked[image] = 1;
-    for (int i = 0; i < (int)m_neighbors[image].size(); ++i)
-      blocked[m_neighbors[image][i]] = 1;
-
-    m_timages[cluster].push_back(image);
-    // m_score, m_cluster, m_satisfied, m_lacks are updated in
-    // setClusters.  So, no need to update anything.
-  }
-  
-  for (int i = 0; i < (int)m_timages.size(); ++i)
-    sort(m_timages[i].begin(), m_timages[i].end());
-  
-  return accumulate(addnum.begin(), addnum.end(), 0);
-}
-
-int Cbundle::totalNum(void) const{
-  int totalnum = 0;
-  for (int c = 0; c < (int)m_timages.size(); ++c)
-    totalnum += (int)m_timages[c].size();
-  return totalnum;
-}
-
-int Cbundle::my_isIntersect(const std::vector<int>& lhs,
-                            const std::vector<int>& rhs) {
-  vector<int>::const_iterator b0 = lhs.begin();
-  vector<int>::const_iterator e0 = lhs.end();
-
-  vector<int>::const_iterator b1 = rhs.begin();
-  vector<int>::const_iterator e1 = rhs.end();
-
-  while (1) {
-    if (b0 == e0)
-      return 0;  
-    if (b1 == e1)
-      return 0;
-
-    if (*b0 == *b1) {
-      return 1;
-    }
-    else if (*b0 < *b1)
-      ++b0;
-    else
-      ++b1;
-  }
-}
-
-void Cbundle::mymerge(const std::vector<int>& lhs,
-                      const std::vector<int>& rhs,
-                      std::vector<int>& output) {
-  vector<int>::const_iterator b0 = lhs.begin();
-  vector<int>::const_iterator e0 = lhs.end();
-
-  vector<int>::const_iterator b1 = rhs.begin();
-  vector<int>::const_iterator e1 = rhs.end();
-
-  while (1) {
-    if (b0 == e0) {
-      output.insert(output.end(), b1, e1);
-      break;
-    }
-    if (b1 == e1) {
-      output.insert(output.end(), b0, e0);
-      break;
-    }
-
-    if (*b0 == *b1) {
-      output.push_back(*b0);
-      ++b0;
-      ++b1;
-    }
-    else if (*b0 < *b1) {
-      output.push_back(*b0);
-      ++b0;
-    }
-    else {
-      output.push_back(*b1);
-      ++b1;
-    }
-  }
-}
-
-// Calculates log2 of number.
-template <typename T>
-T Log2( T n )
-{
-    // log(n)/log(2) is log2.
-    return log( n ) / log( T(2) );
-}
-
-void Cbundle::setWidthsHeightsLevels(void) {
-  m_widths.resize(m_cnum);
-  m_heights.resize(m_cnum);
-  m_levels.resize(m_cnum);
-
-  // Determine level. SfM was done on 2M pixels.
-  for (int c = 0; c < m_cnum; ++c) {
-    m_levels[c] = m_dlevel;
-
-    m_widths[c] = m_pss.m_photos[c].getWidth(m_levels[c]);
-    m_heights[c] = m_pss.m_photos[c].getHeight(m_levels[c]);
-  }
-}
-
- 
-float Cbundle::computeScore2(const Vec4f& coord,
-                             const std::vector<int>& images) const {
-  vector<int> uimages;
-  return computeScore2(coord, images, uimages);
-}
-
-float Cbundle::computeScore2(const Vec4f& coord,
-                             const std::vector<int>& images,
-                             std::vector<int>& uimages) const {
-  const int inum = (int)images.size();
-  uimages.clear();
-  if (inum < 2)
-    return -1.0f;
-  
-  vector<Vec4f> rays;    rays.resize(inum);
-  vector<float> scales;  scales.resize(inum);
-  
-  for (int r = 0; r < inum; ++r) {
-    rays[r] = m_pss.m_photos[images[r]].m_center - coord;
-    unitize(rays[r]);
-    
-    scales[r] = 1.0f / m_pss.m_photos[images[r]].getScale(coord, 0);
-  }
-
-  // Find the best pair of images
-  Vec2i bestpair;
-  float bestscore = -INT_MAX/2;
-  for (int i = 0; i < inum; ++i)
-    for (int j = i+1; j < inum; ++j) {
-      const float ftmp =
-        angleScore(rays[i], rays[j]) * scales[i] * scales[j];
-      if (bestscore < ftmp) {
-        bestscore = ftmp;
-        bestpair = Vec2i(i, j);
-      }
-    }
-
-  vector<int> inout;  inout.resize(inum, 1);
-  inout[bestpair[0]] = 0;         inout[bestpair[1]] = 0;
-  uimages.push_back(bestpair[0]);   uimages.push_back(bestpair[1]);
-
-  for (int t = 2; t < min(m_tau, inum); ++t) {
-    // Add the best new image
-    int ansid = -1;
-    float ansscore = -INT_MAX/2;
-    for (int j = 0; j < inum; ++j) {
-      if (inout[j] == 0)
-        continue;
-
-      float score = 0.0f;
-      for (int k = 0; k < (int)uimages.size(); ++k) {
-        const int iid = uimages[k];
-        score += angleScore(rays[j], rays[iid]) * scales[j] * scales[iid];
-      }
-      if (ansscore < score) {
-        ansscore = score;
-        ansid = j;
-      }
-    }
-    
-    if (ansid == -1) {
-      cerr << "Impossible 2 in compureScore" << endl;      exit (1);
-    }
-    
-    inout[ansid] = 0;
-    uimages.push_back(ansid);
-    bestscore += ansscore;
-  }
-
-  for (int i = 0; i < (int)uimages.size(); ++i)
-    uimages[i] = images[uimages[i]];
-  
-  return bestscore;
-}
-
-float Cbundle::computeScore2(const Vec4f& coord,
-                             const std::vector<int>& images,
-                             const int index) const {
-  vector<int> uimages;
-  const int inum = (int)images.size();
-  if (inum < 2)
-    return -1.0f;
-  
-  vector<Vec4f> rays;    rays.resize(inum);
-  vector<float> scales;  scales.resize(inum);
-  
-  for (int r = 0; r < inum; ++r) {
-    rays[r] = m_pss.m_photos[images[r]].m_center - coord;
-    unitize(rays[r]);
-    
-    scales[r] = 1.0f / m_pss.m_photos[images[r]].getScale(coord, 0);
-  }
-
-  float bestscore = 0.0f;
-  vector<int> inout;  inout.resize(inum, 1);
-  inout[index] = 0;
-  uimages.push_back(index);
-
-  for (int t = 1; t < min(m_tau, inum); ++t) {
-    // Add the best new image
-    int ansid = -1;
-    float ansscore = -INT_MAX/2;
-    for (int j = 0; j < inum; ++j) {
-      if (inout[j] == 0)
-        continue;
-
-      float score = 0.0f;
-      for (int k = 0; k < (int)uimages.size(); ++k) {
-        const int iid = uimages[k];
-        score += angleScore(rays[j], rays[iid]) * scales[j] * scales[iid];
-      }
-      if (ansscore < score) {
-        ansscore = score;
-        ansid = j;
-      }
-    }
-    
-    if (ansid == -1) {
-      cerr << "Impossible 2 in compureScore" << endl;      exit (1);
-    }
-    
-    inout[ansid] = 0;
-    uimages.push_back(ansid);
-    bestscore += ansscore;
-  }
-  return bestscore;
-}
-
-void Cbundle::writeVis(void) {
-  ofstream ofstr;
-  char buffer[1024];
-  sprintf(buffer, "%svis.dat", m_prefix.c_str());
-  
-  ofstr.open(buffer);
-  ofstr << "VISDATA" << endl;
-  ofstr << m_cnum << endl;
-
-  int numer = 0;
-  int denom = 0;
-  
-  for (int c = 0; c < m_cnum; ++c) {
-    if (m_removed[c])
-      ofstr << c << ' ' << 0 << endl;
-    else {
-      ofstr << c << ' ' << (int)m_neighbors[c].size() << "  ";
-      for (int i = 0; i < (int)m_neighbors[c].size(); ++i)
-        ofstr << m_neighbors[c][i] << ' ';
-      ofstr << endl;
-
-      numer += (int)m_neighbors[c].size();
-      ++denom;
-    }
-  }
-  ofstr.close();
-
-
-  cout << numer / (float)denom << " images in vis on the average" << endl;
-  
-}
-
-void Cbundle::writeCameraCenters(void) {
-  for (int i = 0; i < (int)m_timages.size(); ++i) {
-    char buffer[1024];
-    sprintf(buffer, "%scenters-%04d.ply", m_prefix.c_str(), i);
-    
-    ofstream ofstr;
-    ofstr.open(buffer);
-    ofstr << "ply" << endl
-          << "format ascii 1.0" << endl
-          << "element vertex " << (int)m_timages[i].size() << endl
-          << "property float x" << endl
-          << "property float y" << endl
-          << "property float z" << endl
-          << "end_header" << endl;
-    for (int j = 0; j < (int)m_timages[i].size(); ++j) {
-      const int c = m_timages[i][j];
-      ofstr << m_pss.m_photos[c].m_center[0] << ' '
-            << m_pss.m_photos[c].m_center[1] << ' '
-            << m_pss.m_photos[c].m_center[2] << endl;
-    }
-    ofstr.close();
-  }
-
-  {
-    char buffer[1024];
-    sprintf(buffer, "%scenters-all.ply", m_prefix.c_str());
-    ofstream ofstr;
-    ofstr.open(buffer);
-    
-    ofstr << "ply" << endl
-          << "format ascii 1.0" << endl
-          << "element vertex " << m_cnum << endl
-          << "property float x" << endl
-          << "property float y" << endl
-          << "property float z" << endl
-          << "property uchar diffuse_red" << endl
-          << "property uchar diffuse_green" << endl
-          << "property uchar diffuse_blue" << endl
-          << "end_header" << endl;
-    for (int c = 0; c < m_cnum; ++c) {
-      ofstr << m_pss.m_photos[c].m_center[0] << ' '
-            << m_pss.m_photos[c].m_center[1] << ' '
-            << m_pss.m_photos[c].m_center[2] << ' ';
-      
-      if (m_removed[c])
-        ofstr << "0 255 0" << endl;
-      else
-        ofstr << "255 0 255" << endl;
-    }
-    ofstr.close();
-  }
-}
-
-void Cbundle::writeGroups(void) {
-  char buffer[1024];
-  sprintf(buffer, "%sske.dat", m_prefix.c_str());
-  ofstream ofstr;
-  ofstr.open(buffer);
-  ofstr << "SKE" << endl
-        << m_cnum << ' ' << (int)m_timages.size() << endl;
-  
-  for (int c = 0; c < (int)m_timages.size(); ++c) {
-    ofstr << (int)m_timages[c].size() << ' ' << (int)m_oimages[c].size() << endl;
-
-    for (int i = 0; i < (int)m_timages[c].size(); ++i)
-      ofstr << m_timages[c][i] << ' ';
-    ofstr << endl;
-    for (int i = 0; i < (int)m_oimages[c].size(); ++i)
-      ofstr << m_oimages[c][i] << ' ';
-    ofstr << endl;
-  }
-}
-
-void Cbundle::startTimer(void) {
-  time(&m_tv); 
-}
-
-time_t Cbundle::curTimer(void) {
-  time(&m_curtime); 
-  return m_tv - m_curtime;
-}