Rework detector API and implement Harris detector

Switch the detector API to a single function which accepts
a float image and detector options. This makes usage of
feature detection more unified across different algorithms.

Options structure is pretty much straightforward and contains
detector to be used and all the detector-specific settings.

Also implemented Harris feature detection algorithm which
is not as fast as FAST one but is expected to detect more
robust feature points. It is also likely that less features
are detected, but better quality than quantity.

Blender will now use Harris detector by default, later we'll
remove FAST detector.

11 files changed, 514 insertions, 189 deletions

diff --git a/extern/libmv/CMakeLists.txt b/extern/libmv/CMakeLists.txt
index 80a212b64f8..ef7e49b4d81 100644
--- a/extern/libmv/CMakeLists.txt
+++ b/extern/libmv/CMakeLists.txt
@@ -103,6 +103,7 @@ if(WITH_LIBMV)
 		libmv/image/array_nd.h
 		libmv/image/convolve.h
 		libmv/image/correlation.h
+		libmv/image/image_converter.h
 		libmv/image/image.h
 		libmv/image/sample.h
 		libmv/image/tuple.h
diff --git a/extern/libmv/ChangeLog b/extern/libmv/ChangeLog
index c34ea066172..225138b0166 100644
--- a/extern/libmv/ChangeLog
+++ b/extern/libmv/ChangeLog
@@ -1,3 +1,30 @@
+commit b63b8d6989f460fda7d963a2c8b21e8ffa6f8066
+Author: Sergey Sharybin <sergey.vfx@gmail.com>
+Date:   Fri Jan 24 19:32:34 2014 +0600
+
+    Rework detector API and implement Harris detector
+    
+    Switch the detector API to a single function which accepts
+    a float image and detector options. This makes usage of
+    feature detection more unified across different algorithms.
+    
+    Options structure is pretty much straightforward and contains
+    detector to be used and all the detector-specific settings.
+    
+    Also implemented Harris feature detection algorithm which
+    is not as fast as FAST one but is expected to detect more
+    robust feature points.
+    
+    Reviewers: keir
+    
+    Differential Revision: https://developer.blender.org/D258
+
+commit 6458915f64fceba108c5279b7320ca8c76e8a742
+Author: Sergey Sharybin <sergey.vfx@gmail.com>
+Date:   Fri Jan 24 19:14:18 2014 +0600
+
+    Add arcanist configuration file
+
 commit 0a69fadadc5aacbd339f839ac5bd12c3571278b1
 Author: Sergey Sharybin <sergey.vfx@gmail.com>
 Date:   Thu Jan 9 15:50:11 2014 +0600
@@ -659,29 +686,3 @@ Date:   Sat Apr 6 00:38:40 2013 +0600
     
     Also fixed variable shadowing which lead to wrong
     markers positions (were either zero or undefined).
-
-commit a597ca9291ff7c6e209e0500945568afa84d8b4e
-Author: Sergey Sharybin <sergey.vfx@gmail.com>
-Date:   Sat Apr 6 00:12:12 2013 +0600
-
-    Camera intrinsics unit tests fixes
-    
-    - Existing test ApplyIsInvertibleSimple was not
-      doing right thing - distortion model used there
-      was ininvertible.
-    
-      Tweaked parameters in a way model is invertible now.
-    
-    - Added some own tests which tests;
-    
-      * Principal point always maps from pixel space to
-        zero normalized position.
-    
-      * Some basic tests to check whether individual
-        apply/invert works correct.
-
-commit a9a3016c7beb02e62655221d46783b206010d04c
-Author: Sergey Sharybin <sergey.vfx@gmail.com>
-Date:   Thu Apr 4 02:59:58 2013 +0600
-
-    Suppress strict compiler warnings in glags/glog libraries
diff --git a/extern/libmv/files.txt b/extern/libmv/files.txt
index e5297c7cced..a639ea3afa5 100644
--- a/extern/libmv/files.txt
+++ b/extern/libmv/files.txt
@@ -7,6 +7,7 @@ libmv/image/array_nd.h
 libmv/image/convolve.cc
 libmv/image/convolve.h
 libmv/image/correlation.h
+libmv/image/image_converter.h
 libmv/image/image.h
 libmv/image/sample.h
 libmv/image/tuple.h
diff --git a/extern/libmv/libmv-capi.cc b/extern/libmv/libmv-capi.cc
index 474501afbcc..9a536f718aa 100644
--- a/extern/libmv/libmv-capi.cc
+++ b/extern/libmv/libmv-capi.cc
@@ -58,6 +58,14 @@
 #include "libmv/simple_pipeline/reconstruction_scale.h"
 #include "libmv/simple_pipeline/keyframe_selection.h"
 
+#include "libmv/multiview/homography.h"
+
+#include "libmv/image/convolve.h"
+
+#ifdef _MSC_VER
+#  define snprintf _snprintf
+#endif
+
 struct libmv_Reconstruction {
 	libmv::EuclideanReconstruction reconstruction;
 
@@ -69,7 +77,7 @@ struct libmv_Reconstruction {
 };
 
 struct libmv_Features {
-	int count, margin;
+	int count;
 	libmv::Feature *features;
 };
 
@@ -107,6 +115,21 @@ void libmv_setLoggingVerbosity(int verbosity)
 
 /* ************ Utility ************ */
 
+static void byteBufToImage(const unsigned char *buf, int width, int height, int channels, libmv::FloatImage *image)
+{
+	int x, y, k, a = 0;
+
+	image->Resize(height, width, channels);
+
+	for (y = 0; y < height; y++) {
+		for (x = 0; x < width; x++) {
+			for (k = 0; k < channels; k++) {
+				(*image)(y, x, k) = (float)buf[a++] / 255.0f;
+			}
+		}
+	}
+}
+
 static void floatBufToImage(const float *buf, int width, int height, int channels, libmv::FloatImage *image)
 {
 	int x, y, k, a = 0;
@@ -833,63 +856,88 @@ struct libmv_CameraIntrinsics *libmv_reconstructionExtractIntrinsics(struct libm
 
 /* ************ Feature detector ************ */
 
-struct libmv_Features *libmv_detectFeaturesFAST(const unsigned char *data,
-                                                int width, int height, int stride,
-                                                int margin, int min_trackness, int min_distance)
+static libmv_Features *libmv_featuresFromVector(
+	const libmv::vector<libmv::Feature> &features)
 {
-	libmv::Feature *features = NULL;
-	std::vector<libmv::Feature> v;
 	struct libmv_Features *libmv_features = LIBMV_STRUCT_NEW(libmv_Features, 1);
-	int i = 0, count;
-
-	if (margin) {
-		data += margin * stride+margin;
-		width -= 2 * margin;
-		height -= 2 * margin;
-	}
-
-	v = libmv::DetectFAST(data, width, height, stride, min_trackness, min_distance);
-
-	count = v.size();
-
+	int count = features.size();
 	if (count) {
-		features = LIBMV_STRUCT_NEW(libmv::Feature, count);
+		libmv_features->features = LIBMV_STRUCT_NEW(libmv::Feature, count);
 
-		for(std::vector<libmv::Feature>::iterator it = v.begin(); it != v.end(); it++) {
-			features[i++] = *it;
+		for (int i = 0; i < count; i++) {
+			libmv_features->features[i] = features.at(i);
 		}
 	}
+	else {
+		libmv_features->features = NULL;
+	}
 
-	libmv_features->features = features;
 	libmv_features->count = count;
-	libmv_features->margin = margin;
 
-	return (struct libmv_Features *)libmv_features;
+	return libmv_features;
 }
 
-struct libmv_Features *libmv_detectFeaturesMORAVEC(const unsigned char *data,
-                                                   int width, int height, int stride,
-                                                   int margin, int count, int min_distance)
+static void libmv_convertDetectorOptions(libmv_DetectOptions *options,
+                                         libmv::DetectOptions *detector_options)
 {
-	libmv::Feature *features = NULL;
-	struct libmv_Features *libmv_features = LIBMV_STRUCT_NEW(libmv_Features, 1);
+	switch (options->detector) {
+#define LIBMV_CONVERT(the_detector) \
+	case LIBMV_DETECTOR_ ## the_detector: \
+		detector_options->type = libmv::DetectOptions::the_detector; \
+		break;
+		LIBMV_CONVERT(FAST)
+		LIBMV_CONVERT(MORAVEC)
+		LIBMV_CONVERT(HARRIS)
+#undef LIBMV_CONVERT
+	}
+	detector_options->margin = options->margin;
+	detector_options->min_distance = options->min_distance;
+	detector_options->fast_min_trackness = options->fast_min_trackness;
+	detector_options->moravec_max_count = options->moravec_max_count;
+	detector_options->moravec_pattern = options->moravec_pattern;
+	detector_options->harris_threshold = options->harris_threshold;
+}
 
-	if (count) {
-		if (margin) {
-			data += margin * stride+margin;
-			width -= 2 * margin;
-			height -= 2 * margin;
-		}
+struct libmv_Features *libmv_detectFeaturesByte(const unsigned char *image_buffer,
+                                                int width, int height, int channels,
+                                                libmv_DetectOptions *options)
+{
+	// Prepare the image.
+	libmv::FloatImage image;
+	byteBufToImage(image_buffer, width, height, channels, &image);
 
-		features = LIBMV_STRUCT_NEW(libmv::Feature, count);
-		libmv::DetectMORAVEC(data, stride, width, height, features, &count, min_distance, NULL);
-	}
+	// Configure detector.
+	libmv::DetectOptions detector_options;
+	libmv_convertDetectorOptions(options, &detector_options);
 
-	libmv_features->count = count;
-	libmv_features->margin = margin;
-	libmv_features->features = features;
+	// Run the detector.
+	libmv::vector<libmv::Feature> detected_features;
+	libmv::Detect(image, detector_options, &detected_features);
 
-	return libmv_features;
+	// Convert result to C-API.
+	libmv_Features *result = libmv_featuresFromVector(detected_features);
+	return result;
+}
+
+struct libmv_Features *libmv_detectFeaturesFloat(const float *image_buffer,
+                                                 int width, int height, int channels,
+                                                 libmv_DetectOptions *options)
+{
+	// Prepare the image.
+	libmv::FloatImage image;
+	floatBufToImage(image_buffer, width, height, channels, &image);
+
+	// Configure detector.
+	libmv::DetectOptions detector_options;
+	libmv_convertDetectorOptions(options, &detector_options);
+
+	// Run the detector.
+	libmv::vector<libmv::Feature> detected_features;
+	libmv::Detect(image, detector_options, &detected_features);
+
+	// Convert result to C-API.
+	libmv_Features *result = libmv_featuresFromVector(detected_features);
+	return result;
 }
 
 void libmv_featuresDestroy(struct libmv_Features *libmv_features)
@@ -910,8 +958,8 @@ void libmv_getFeature(const struct libmv_Features *libmv_features, int number, d
 {
 	libmv::Feature feature = libmv_features->features[number];
 
-	*x = feature.x + libmv_features->margin;
-	*y = feature.y + libmv_features->margin;
+	*x = feature.x;
+	*y = feature.y;
 	*score = feature.score;
 	*size = feature.size;
 }
diff --git a/extern/libmv/libmv-capi.h b/extern/libmv/libmv-capi.h
index 13cc3ae8499..d183bc4cd41 100644
--- a/extern/libmv/libmv-capi.h
+++ b/extern/libmv/libmv-capi.h
@@ -118,10 +118,29 @@ double libmv_reprojectionError(const struct libmv_Reconstruction *libmv_reconstr
 struct libmv_CameraIntrinsics *libmv_reconstructionExtractIntrinsics(struct libmv_Reconstruction *libmv_Reconstruction);
 
 /* Feature detector */
-struct libmv_Features *libmv_detectFeaturesFAST(const unsigned char *data, int width, int height, int stride,
-                                                int margin, int min_trackness, int min_distance);
-struct libmv_Features *libmv_detectFeaturesMORAVEC(const unsigned char *data, int width, int height, int stride,
-                                                   int margin, int count, int min_distance);
+enum {
+	LIBMV_DETECTOR_FAST,
+	LIBMV_DETECTOR_MORAVEC,
+	LIBMV_DETECTOR_HARRIS,
+};
+
+typedef struct libmv_DetectOptions {
+  int detector;
+  int margin;
+  int min_distance;
+  int fast_min_trackness;
+  int moravec_max_count;
+  unsigned char *moravec_pattern;
+  double harris_threshold;
+} libmv_DetectOptions;
+
+struct libmv_Features *libmv_detectFeaturesByte(const unsigned char *image_buffer,
+                                                int width, int height, int channels,
+                                                libmv_DetectOptions *options);
+struct libmv_Features *libmv_detectFeaturesFloat(const float *image_buffer,
+                                                 int width, int height, int channels,
+                                                 libmv_DetectOptions *options);
+
 void libmv_featuresDestroy(struct libmv_Features *libmv_features);
 int libmv_countFeatures(const struct libmv_Features *libmv_features);
 void libmv_getFeature(const struct libmv_Features *libmv_features, int number, double *x, double *y, double *score,
diff --git a/extern/libmv/libmv-capi_stub.cc b/extern/libmv/libmv-capi_stub.cc
index bda9605b422..bd5d16c9077 100644
--- a/extern/libmv/libmv-capi_stub.cc
+++ b/extern/libmv/libmv-capi_stub.cc
@@ -146,15 +146,16 @@ void libmv_reconstructionDestroy(struct libmv_Reconstruction * /*libmv_reconstru
 
 /* ************ feature detector ************ */
 
-struct libmv_Features *libmv_detectFeaturesFAST(const unsigned char * /*data*/, int /*width*/, int /*height*/,
-                                                int /*stride*/, int /*margin*/, int /*min_trackness*/,
-                                                int /*min_distance*/)
+struct libmv_Features *libmv_detectFeaturesByte(const unsigned char */*image_buffer*/,
+                                                int /*width*/, int /*height*/, int /*channels*/,
+                                                libmv_DetectOptions */*options*/)
 {
 	return NULL;
 }
 
-struct libmv_Features *libmv_detectFeaturesMORAVEC(const unsigned char * /*data*/, int /*width*/, int /*height*/,
-                                                   int /*stride*/, int /*margin*/, int /*count*/, int /*min_distance*/)
+struct libmv_Features *libmv_detectFeaturesFloat(const float */*image_buffer*/,
+                                                 int /*width*/, int /*height*/, int /*channels*/,
+                                                 libmv_DetectOptions */*options*/)
 {
 	return NULL;
 }
diff --git a/extern/libmv/libmv/image/array_nd.h b/extern/libmv/libmv/image/array_nd.h
index f40fadcbde7..c5099f24d7b 100644
--- a/extern/libmv/libmv/image/array_nd.h
+++ b/extern/libmv/libmv/image/array_nd.h
@@ -471,6 +471,23 @@ void MultiplyElements(const ArrayND<TA, 3> &a,
   }
 }
 
+template <typename TA, typename TB, typename TC>
+void MultiplyElements(const Array3D<TA> &a,
+                      const Array3D<TB> &b,
+                      Array3D<TC> *c) {
+  // Specialization for N==3
+  c->ResizeLike(a);
+  assert(a.Shape(0) == b.Shape(0));
+  assert(a.Shape(1) == b.Shape(1));
+  assert(a.Shape(2) == b.Shape(2));
+  for (int i = 0; i < a.Shape(0); ++i) {
+    for (int j = 0; j < a.Shape(1); ++j) {
+      for (int k = 0; k < a.Shape(2); ++k) {
+        (*c)(i, j, k) = TC(a(i, j, k) * b(i, j, k));
+      }
+    }
+  }
+}
 
 }  // namespace libmv
 
diff --git a/extern/libmv/libmv/image/convolve.cc b/extern/libmv/libmv/image/convolve.cc
index 4682c619f04..464043581d2 100644
--- a/extern/libmv/libmv/image/convolve.cc
+++ b/extern/libmv/libmv/image/convolve.cc
@@ -182,6 +182,23 @@ void ConvolveGaussian(const Array3Df &in,
   ConvolveHorizontal(tmp, kernel, out_pointer);
 }
 
+void ImageDerivatives(const Array3Df &in,
+                      double sigma,
+                      Array3Df *gradient_x,
+                      Array3Df *gradient_y) {
+  Vec kernel, derivative;
+  ComputeGaussianKernel(sigma, &kernel, &derivative);
+  Array3Df tmp;
+
+  // Compute first derivative in x.
+  ConvolveVertical(in, kernel, &tmp);
+  ConvolveHorizontal(tmp, derivative, gradient_x);
+
+  // Compute first derivative in y.
+  ConvolveHorizontal(in, kernel, &tmp);
+  ConvolveVertical(tmp, derivative, gradient_y);
+}
+
 void BlurredImageAndDerivatives(const Array3Df &in,
                                 double sigma,
                                 Array3Df *blurred_image,
diff --git a/extern/libmv/libmv/image/image_converter.h b/extern/libmv/libmv/image/image_converter.h
new file mode 100644
index 00000000000..b3a3fa2bf8c
--- /dev/null
+++ b/extern/libmv/libmv/image/image_converter.h
@@ -0,0 +1,81 @@
+// Copyright (c) 2009 libmv authors.
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to
+// deal in the Software without restriction, including without limitation the
+// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
+// sell copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in
+// all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
+// IN THE SOFTWARE.
+
+#ifndef LIBMV_IMAGE_IMAGE_CONVERTER_H
+#define LIBMV_IMAGE_IMAGE_CONVERTER_H
+
+#include "libmv/image/array_nd.h"
+
+namespace libmv {
+
+// The factor comes from http://www.easyrgb.com/
+// RGB to XYZ : Y is the luminance channel
+// var_R * 0.2126 + var_G * 0.7152 + var_B * 0.0722
+template<typename T>
+inline T RGB2GRAY(const T r, const T g, const T b) {
+  return static_cast<T>(r * 0.2126 + g * 0.7152 + b * 0.0722);
+}
+
+/*
+// Specialization for the uchar type. (that do not want to work...)
+template<>
+inline unsigned char RGB2GRAY<unsigned char>(const unsigned char r,
+                                             const unsigned char g,
+                                             const unsigned char b) {
+  return (unsigned char)(r * 0.2126 + g * 0.7152 + b * 0.0722 +0.5);
+}*/
+
+template<class ImageIn, class ImageOut>
+void Rgb2Gray(const ImageIn &imaIn, ImageOut *imaOut) {
+  // It is all fine to cnvert RGBA image here as well,
+  // all the additional channels will be nicely ignored.
+  assert(imaIn.Depth() >= 3);
+
+  imaOut->Resize(imaIn.Height(), imaIn.Width(), 1);
+  // Convert each RGB pixel into Gray value (luminance)
+
+  for (int j = 0; j < imaIn.Height(); ++j) {
+    for (int i = 0; i < imaIn.Width(); ++i)  {
+      (*imaOut)(j, i) = RGB2GRAY(imaIn(j, i, 0) , imaIn(j, i, 1), imaIn(j, i, 2));
+    }
+  }
+}
+
+// Convert given float image to an unsigned char array of pixels.
+template<class Image>
+unsigned char *FloatImageToUCharArray(const Image &image) {
+  unsigned char *buffer =
+      new unsigned char[image.Width() * image.Height() * image.Depth()];
+
+  for (int y = 0; y < image.Height(); y++) {
+    for (int x = 0; x < image.Width(); x++)  {
+      for (int d = 0; d < image.Depth(); d++)  {
+        int index = (y * image.Width() + x) * image.Depth() + d;
+        buffer[index] = 255.0 * image(y, x, d);
+      }
+    }
+  }
+
+  return buffer;
+}
+
+}  // namespace libmv
+
+#endif  // LIBMV_IMAGE_IMAGE_CONVERTER_H
diff --git a/extern/libmv/libmv/simple_pipeline/detect.cc b/extern/libmv/libmv/simple_pipeline/detect.cc
index 27639e958a0..1ddb23dcd28 100644
--- a/extern/libmv/libmv/simple_pipeline/detect.cc
+++ b/extern/libmv/libmv/simple_pipeline/detect.cc
@@ -22,96 +22,133 @@
 **
 ****************************************************************************/
 
-#include "libmv/simple_pipeline/detect.h"
-#include <third_party/fast/fast.h>
 #include <stdlib.h>
 #include <memory.h>
+#include <queue>
+#include <third_party/fast/fast.h>
+
+#include "libmv/base/scoped_ptr.h"
+#include "libmv/image/array_nd.h"
+#include "libmv/image/image_converter.h"
+#include "libmv/image/convolve.h"
+#include "libmv/logging/logging.h"
+#include "libmv/simple_pipeline/detect.h"
 
 #ifdef __SSE2__
-#include <emmintrin.h>
+#  include <emmintrin.h>
 #endif
 
 namespace libmv {
 
-typedef unsigned int uint;
+namespace {
+
+class FeatureComparison {
+ public:
+  bool operator() (const Feature &left, const Feature &right) const {
+    return right.score > left.score;
+  }
+};
+
+// Filter the features so there are no features closer than
+// minimal distance to each other.
+// This is a naive implementation with O(n^2) asymptotic.
+void FilterFeaturesByDistance(const vector<Feature> &all_features,
+                              int min_distance,
+                              vector<Feature> *detected_features) {
+  const int min_distance_squared = min_distance * min_distance;
+
+  // Use priority queue to sort the features by their score.
+  //
+  // Do this on copy of the input features to prevent possible
+  // distortion in callee function behavior.
+  std::priority_queue<Feature,
+                      std::vector<Feature>,
+                      FeatureComparison> priority_features;
+
+  for (int i = 0; i < all_features.size(); i++) {
+    priority_features.push(all_features.at(i));
+  }
+
+  while (!priority_features.empty()) {
+    bool ok = true;
+    Feature a = priority_features.top();
+
+    for (int i = 0; i < detected_features->size(); i++) {
+      Feature &b = detected_features->at(i);
+      if (Square(a.x - b.x) + Square(a.y - b.y) < min_distance_squared) {
+        ok = false;
+        break;
+      }
+    }
 
-static int featurecmp(const void *a_v, const void *b_v) {
-  Feature *a = (Feature*)a_v;
-  Feature *b = (Feature*)b_v;
+    if (ok) {
+      detected_features->push_back(a);
+    }
 
-  return b->score - a->score;
+    priority_features.pop();
+  }
 }
 
-std::vector<Feature> DetectFAST(const unsigned char* data,
-                                int width, int height,
-                                int stride,
-                                int min_trackness,
-                                int min_distance) {
-  std::vector<Feature> features;
+void DetectFAST(const FloatImage &grayscale_image,
+                const DetectOptions &options,
+                vector<Feature> *detected_features) {
+  const int min_distance = options.min_distance;
+  const int min_trackness = options.fast_min_trackness;
+  const int margin = options.margin;
+  const int width = grayscale_image.Width() - 2 * margin;
+  const int height = grayscale_image.Width() - 2 * margin;
+  const int stride = grayscale_image.Width();
+
+  // TODO(sergey): Use scoped_array to guard de-allocation of the array.
+  // Same goes to `scores` and `all` arrays here.
+  unsigned char *byte_image = FloatImageToUCharArray(grayscale_image);
+  const int byte_image_offset = margin * stride + margin;
+
   // TODO(MatthiasF): Support targetting a feature count (binary search trackness)
   int num_features;
-  xy* all = fast9_detect(data, width, height,
-                         stride, min_trackness, &num_features);
+  xy *all = fast9_detect(byte_image + byte_image_offset,
+                         width,
+                         height,
+                         stride,
+                         min_trackness,
+                         &num_features);
   if (num_features == 0) {
     free(all);
-    return features;
+    delete [] byte_image;
+    return;
   }
-  int* scores = fast9_score(data, stride, all, num_features, min_trackness);
+  int *scores = fast9_score(byte_image + byte_image_offset,
+                            stride,
+                            all,
+                            num_features,
+                            min_trackness);
   // TODO(MatthiasF): merge with close feature suppression
-  xy* nonmax = nonmax_suppression(all, scores, num_features, &num_features);
+  xy *nonmax = nonmax_suppression(all, scores, num_features, &num_features);
   free(all);
+  delete [] byte_image;
   // Remove too close features
   // TODO(MatthiasF): A resolution independent parameter would be better than
   // distance e.g. a coefficient going from 0 (no minimal distance) to 1
   // (optimal circle packing)
   // FIXME(MatthiasF): this method will not necessarily give all maximum markers
   if (num_features) {
-    Feature *all_features = new Feature[num_features];
-
+    vector<Feature> all_features;
     for (int i = 0; i < num_features; ++i) {
-      Feature a = { (float)nonmax[i].x, (float)nonmax[i].y, (float)scores[i], 0 };
-      all_features[i] = a;
-    }
-
-    qsort((void *)all_features, num_features, sizeof(Feature), featurecmp);
-
-    features.reserve(num_features);
-
-    int prev_score = all_features[0].score;
-    const int min_distance_squared = min_distance * min_distance;
-    for (int i = 0; i < num_features; ++i) {
-      bool ok = true;
-      Feature a = all_features[i];
-      if (a.score>prev_score)
-        abort();
-      prev_score = a.score;
-
-      // compare each feature against filtered set
-      for (int j = 0; j < features.size(); j++) {
-        Feature& b = features[j];
-        if ((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y) < min_distance_squared) {
-          // already a nearby feature
-          ok = false;
-          break;
-        }
-      }
-
-      if (ok) {
-        // add the new feature
-        features.push_back(a);
-      }
+      Feature new_feature = {(float)nonmax[i].x + margin,
+                             (float)nonmax[i].y + margin,
+                             (float)scores[i],
+                             0};
+      all_features.push_back(new_feature);
     }
-
-    delete [] all_features;
+    FilterFeaturesByDistance(all_features, min_distance, detected_features);
   }
   free(scores);
   free(nonmax);
-  return features;
 }
 
 #ifdef __SSE2__
-static uint SAD(const ubyte* imageA, const ubyte* imageB,
-                int strideA, int strideB) {
+static unsigned int SAD(const ubyte* imageA, const ubyte* imageB,
+                        int strideA, int strideB) {
   __m128i a = _mm_setzero_si128();
   for (int i = 0; i < 16; i++) {
     a = _mm_adds_epu16(a,
@@ -121,9 +158,9 @@ static uint SAD(const ubyte* imageA, const ubyte* imageB,
   return _mm_extract_epi16(a, 0) + _mm_extract_epi16(a, 4);
 }
 #else
-static uint SAD(const ubyte* imageA, const ubyte* imageB,
-                int strideA, int strideB) {
-  uint sad = 0;
+static unsigned int SAD(const ubyte* imageA, const ubyte* imageB,
+                        int strideA, int strideB) {
+  unsigned int sad = 0;
   for (int i = 0; i < 16; i++) {
     for (int j = 0; j < 16; j++) {
       sad += abs((int)imageA[i*strideA+j] - imageB[i*strideB+j]);
@@ -133,11 +170,20 @@ static uint SAD(const ubyte* imageA, const ubyte* imageB,
 }
 #endif
 
-void DetectMORAVEC(const ubyte* image,
-                   int stride, int width, int height,
-                   Feature* detected, int* count,
-                   int distance,
-                   ubyte* pattern) {
+void DetectMORAVEC(const FloatImage &grayscale_image,
+                   const DetectOptions &options,
+                   vector<Feature> *detected_features) {
+  const int distance = options.min_distance;
+  const int margin = options.margin;
+  const unsigned char *pattern = options.moravec_pattern;
+  const int count = options.moravec_max_count;
+  const int width = grayscale_image.Width() - 2 * margin;
+  const int height = grayscale_image.Width() - 2 * margin;
+  const int stride = grayscale_image.Width();
+
+  // TODO(sergey): Use scoped_array to guard de-allocation of the array.
+  unsigned char *byte_image = FloatImageToUCharArray(grayscale_image);
+
   unsigned short histogram[256];
   memset(histogram, 0, sizeof(histogram));
   ubyte* scores = new ubyte[width*height];
@@ -145,7 +191,7 @@ void DetectMORAVEC(const ubyte* image,
   const int r = 1;  // radius for self similarity comparison
   for (int y = distance; y < height-distance; y++) {
     for (int x = distance; x < width-distance; x++) {
-      const ubyte* s = &image[y*stride+x];
+      const ubyte* s = &byte_image[y*stride+x];
       int score =  // low self-similarity with overlapping patterns
                    // OPTI: load pattern once
           SAD(s, s-r*stride-r, stride, stride)+SAD(s, s-r*stride, stride, stride)+SAD(s, s-r*stride+r, stride, stride)+
@@ -182,19 +228,105 @@ void DetectMORAVEC(const ubyte* image,
   int min = 255, total = 0;
   for (; min > 0; min--) {
     int h = histogram[min];
-    if (total+h > *count) break;
+    if (total+h > count) break;
     total += h;
   }
-  int i = 0;
   for (int y = 16; y < height-16; y++) {
     for (int x = 16; x < width-16; x++) {
       int s = scores[y*width+x];
       Feature f = { (float)x+8.0f, (float)y+8.0f, (float)s, 16 };
-      if (s > min) detected[i++] = f;
+      if (s > min) {
+        detected_features->push_back(f);
+      }
     }
   }
-  *count = i;
   delete[] scores;
+  delete[] byte_image;
+}
+
+void DetectHarris(const FloatImage &grayscale_image,
+                  const DetectOptions &options,
+                  vector<Feature> *detected_features) {
+  const double alpha = 0.06;
+  const double sigma = 0.9;
+
+  const int min_distance = options.min_distance;
+  const int margin = options.margin;
+  const double threshold = options.harris_threshold;
+
+  FloatImage gradient_x, gradient_y;
+  ImageDerivatives(grayscale_image, sigma, &gradient_x, &gradient_y);
+
+  FloatImage gradient_xx, gradient_yy, gradient_xy;
+  MultiplyElements(gradient_x, gradient_x, &gradient_xx);
+  MultiplyElements(gradient_y, gradient_y, &gradient_yy);
+  MultiplyElements(gradient_x, gradient_y, &gradient_xy);
+
+  FloatImage gradient_xx_blurred,
+             gradient_yy_blurred,
+             gradient_xy_blurred;
+  ConvolveGaussian(gradient_xx, sigma, &gradient_xx_blurred);
+  ConvolveGaussian(gradient_yy, sigma, &gradient_yy_blurred);
+  ConvolveGaussian(gradient_xy, sigma, &gradient_xy_blurred);
+
+  vector<Feature> all_features;
+  for (int y = margin; y < gradient_xx_blurred.Height() - margin; ++y) {
+    for (int x = margin; x < gradient_xx_blurred.Width() - margin; ++x) {
+      // Construct matrix
+      //
+      //  A = [ Ix^2  Ix*Iy ]
+      //      [ Ix*Iy Iy^2  ]
+      Mat2 A;
+      A(0, 0) = gradient_xx_blurred(y, x);
+      A(1, 1) = gradient_yy_blurred(y, x);
+      A(0, 1) = A(1, 0) = gradient_xy_blurred(y, x);
+
+      double detA = A.determinant();
+      double traceA = A.trace();
+      double harris_function = detA - alpha * traceA * traceA;
+      if (harris_function > threshold) {
+        Feature new_feature = {(float)x, (float)y, (float)harris_function, 0.0f};
+        all_features.push_back(new_feature);
+      }
+    }
+  }
+
+  FilterFeaturesByDistance(all_features, min_distance, detected_features);
+}
+
+}  // namespace
+
+DetectOptions::DetectOptions()
+  : type(DetectOptions::HARRIS),
+    margin(0),
+    min_distance(120),
+    fast_min_trackness(128),
+    moravec_max_count(0),
+    moravec_pattern(NULL),
+    harris_threshold(0.0) {}
+
+void Detect(const FloatImage &image,
+            const DetectOptions &options,
+            vector<Feature> *detected_features) {
+  // Currently all the detectors requires image to be grayscale.
+  // Do it here to avoid code duplication.
+  FloatImage grayscale_image;
+  if (image.Depth() != 1) {
+    Rgb2Gray(image, &grayscale_image);
+  } else {
+    // TODO(sergey): Find a way to avoid such image duplication/
+    grayscale_image = image;
+  }
+
+  if (options.type == DetectOptions::FAST) {
+    DetectFAST(grayscale_image, options, detected_features);
+  } else if (options.type == DetectOptions::MORAVEC) {
+    DetectMORAVEC(grayscale_image, options, detected_features);
+  } else if (options.type == DetectOptions::HARRIS) {
+    DetectHarris(grayscale_image, options, detected_features);
+  } else {
+    LOG(FATAL) << "Unknown detector has been passed to featur detection";
+  }
 }
 
 }  // namespace libmv
diff --git a/extern/libmv/libmv/simple_pipeline/detect.h b/extern/libmv/libmv/simple_pipeline/detect.h
index f5bcd07afd7..2c782294d1b 100644
--- a/extern/libmv/libmv/simple_pipeline/detect.h
+++ b/extern/libmv/libmv/simple_pipeline/detect.h
@@ -27,6 +27,9 @@
 
 #include <vector>
 
+#include "libmv/base/vector.h"
+#include "libmv/image/image.h"
+
 namespace libmv {
 
 typedef unsigned char ubyte;
@@ -50,47 +53,51 @@ struct Feature {
   float size;
 };
 
-/*!
-    Detect features in an image.
-
-    You need to input a single channel 8-bit image using pointer to image \a data,
-    \a width, \a height and \a stride (i.e bytes per line).
-
-    You can tweak the count of detected features using \a min_trackness, which is
-    the minimum score to add a feature, and \a min_distance which is the minimal
-    distance accepted between two featuress.
-
-    \note You can binary search over \a min_trackness to get a given feature count.
-
-    \note a way to get an uniform distribution of a given feature count is:
-          \a min_distance = \a width * \a height / desired_feature_count ^ 2
-
-    \return All detected feartures matching given parameters
-*/
-std::vector<Feature> DetectFAST(const unsigned char* data, int width, int height,
-                           int stride, int min_trackness = 128,
-                           int min_distance = 120);
-
-/*!
-    Detect features in an image.
-
-    \a image is a single channel 8-bit image of size \a width x \a height
-
-    \a detected is an array with space to hold \a *count features.
-    \a *count is the maximum count to detect on input and the actual
-    detected count on output.
-
-    \a distance is the minimal distance between detected features.
-
-    if \a pattern is null all good features will be found.
-    if \a pattern is not null only features similar to \a pattern will be found.
+struct DetectOptions {
+  DetectOptions();
+
+  // TODO(sergey): Write descriptions to each of algorithms.
+  enum DetectorType {
+    FAST,
+    MORAVEC,
+    HARRIS,
+  };
+  DetectorType type;
+
+  // Margin in pixels from the image boundary.
+  // No features will be detected within the margin.
+  int margin;
+
+  // Minimal distance between detected features.
+  int min_distance;
+
+  // Minimum score to add a feature. Only used by FAST detector.
+  //
+  // TODO(sergey): Write more detailed documentation about which
+  // units this value is measured in and so on.
+  int fast_min_trackness;
+
+  // Maximum count to detect. Only used by MORAVEC detector.
+  int moravec_max_count;
+
+  // Find only features similar to this pattern. Only used by MORAVEC detector.
+  //
+  // This is an image patch denoted in byte array with dimensions of 16px by 16px
+  // used to filter features by similarity to this patch.
+  unsigned char *moravec_pattern;
+
+  // Threshold value of the Harris function to add new featrue
+  // to the result.
+  double harris_threshold;
+};
 
-    \note \a You can crop the image (to avoid detecting markers near the borders) without copying:
-             image += marginY*stride+marginX, width -= 2*marginX, height -= 2*marginY;
-*/
-void DetectMORAVEC(const ubyte* image, int stride, int width, int height,
-                   Feature* detected, int* count, int distance /*=32*/,
-                   ubyte* pattern /*=0*/);
+// Detect features on a given image using given detector options.
+//
+// Image could have 1-4 channels, it'll be converted to a grayscale
+// by the detector function if needed.
+void Detect(const FloatImage &image,
+            const DetectOptions &options,
+            vector<Feature> *detected_features);
 
 }  // namespace libmv