arrayfire · pavanky · Nov 11, 2015 · Nov 3, 2015 · Nov 3, 2015 · Nov 3, 2015
diff --git a/docs/details/vision.dox b/docs/details/vision.dox
@@ -170,6 +170,30 @@ Template matching is an image processing technique to find small patches of an i
 match a given template image. A more in depth discussion on the topic can be found
 [here](http://en.wikipedia.org/wiki/Template_matching).
 
+=======================================================================
+
+\defgroup cv_func_homography homography
+\ingroup homography_mat
+
+\brief Homography Estimation
+
+Homography estimation find a perspective transform between two sets of 2D points.
+Currently, two methods are supported for the estimation, RANSAC (RANdom SAmple Consensus)
+and LMedS (Least Median of Squares). Both methods work by randomly selecting a subset
+of 4 points of the set of source points, computing the eigenvectors of that set and
+finding the perspective transform. The process is repeated several times, a maximum of
+times given by the value passed to the iterations arguments for RANSAC (for the CPU
+backend, usually less than that, depending on the quality of the dataset, but for CUDA
+and OpenCL backends the transformation will be computed exactly the amount of times
+passed via the iterations parameter), the returned value is the one that matches the
+best number of inliers, which are all of the points that fall within a maximum L2
+distance from the value passed to the inlier_thr argument. For the LMedS case, the
+number of iterations is currently hardcoded to meet the following equation:
+
+\f$ m = \frac{log(1 - P)}{log[1 - {(1 - \epsilon)}^{p}]}\f$,
+
+where \f$ P = 0.99\f$, \f$ \epsilon = 40\%\f$ and \f$ p = 4\f$.
+
 
 
 @}

diff --git a/include/af/defines.h b/include/af/defines.h
@@ -325,6 +325,11 @@ typedef enum {
     AF_FIF_RAW          = 34    ///< FreeImage Enum for RAW Camera Image File
 } af_image_format;
 
+typedef enum {
+    AF_HOMOGRAPHY_RANSAC = 0,   ///< Computes homography using RANSAC
+    AF_HOMOGRAPHY_LMEDS  = 1    ///< Computes homography using Least Median of Squares
+} af_homography_type;
+
 // These enums should be 2^x
 typedef enum {
     AF_BACKEND_DEFAULT = 0,  ///< Default backend order: OpenCL -> CUDA -> CPU

diff --git a/include/af/vision.h b/include/af/vision.h
@@ -39,7 +39,9 @@ class array;
 
     \ingroup cv_func_fast
  */
-AFAPI features fast(const array& in, const float thr=20.0f, const unsigned arc_length=9, const bool non_max=true, const float feature_ratio=0.05, const unsigned edge=3);
+AFAPI features fast(const array& in, const float thr=20.0f, const unsigned arc_length=9,
+                    const bool non_max=true, const float feature_ratio=0.05,
+                    const unsigned edge=3);
 
 #if AF_API_VERSION >= 31
 /**
@@ -68,7 +70,9 @@ AFAPI features fast(const array& in, const float thr=20.0f, const unsigned arc_l
 
     \ingroup cv_func_harris
  */
-AFAPI features harris(const array& in, const unsigned max_corners=500, const float min_response=1e5f, const float sigma=1.f, const unsigned block_size=0, const float k_thr=0.04f);
+AFAPI features harris(const array& in, const unsigned max_corners=500,
+                      const float min_response=1e5f, const float sigma=1.f,
+                      const unsigned block_size=0, const float k_thr=0.04f);
 #endif
 
 /**
@@ -93,7 +97,10 @@ AFAPI features harris(const array& in, const unsigned max_corners=500, const flo
 
     \ingroup cv_func_orb
  */
-AFAPI void orb(features& feat, array& desc, const array& image, const float fast_thr=20.f, const unsigned max_feat=400, const float scl_fctr=1.5f, const unsigned levels=4, const bool blur_img=false);
+AFAPI void orb(features& feat, array& desc, const array& image,
+               const float fast_thr=20.f, const unsigned max_feat=400,
+               const float scl_fctr=1.5f, const unsigned levels=4,
+               const bool blur_img=false);
 
 #if AF_API_VERSION >= 31
 /**
@@ -127,7 +134,10 @@ AFAPI void orb(features& feat, array& desc, const array& image, const float fast
 
     \ingroup cv_func_sift
  */
-AFAPI void sift(features& feat, array& desc, const array& in, const unsigned n_layers=3, const float contrast_thr=0.04f, const float edge_thr=10.f, const float init_sigma=1.6f, const bool double_input=true, const float intensity_scale=0.00390625f, const float feature_ratio=0.05f);
+AFAPI void sift(features& feat, array& desc, const array& in, const unsigned n_layers=3,
+                const float contrast_thr=0.04f, const float edge_thr=10.f,
+                const float init_sigma=1.6f, const bool double_input=true,
+                const float intensity_scale=0.00390625f, const float feature_ratio=0.05f);
 #endif
 
 #if AF_API_VERSION >= 32
@@ -162,7 +172,10 @@ AFAPI void sift(features& feat, array& desc, const array& in, const unsigned n_l
 
     \ingroup cv_func_sift
  */
-AFAPI void gloh(features& feat, array& desc, const array& in, const unsigned n_layers=3, const float contrast_thr=0.04f, const float edge_thr=10.f, const float init_sigma=1.6f, const bool double_input=true, const float intensity_scale=0.00390625f, const float feature_ratio=0.05f);
+AFAPI void gloh(features& feat, array& desc, const array& in, const unsigned n_layers=3,
+                const float contrast_thr=0.04f, const float edge_thr=10.f,
+                const float init_sigma=1.6f, const bool double_input=true,
+                const float intensity_scale=0.00390625f, const float feature_ratio=0.05f);
 #endif
 
 /**
@@ -280,6 +293,37 @@ AFAPI features susan(const array& in,
 AFAPI array dog(const array& in, const int radius1, const int radius2);
 #endif
 
+#if AF_API_VERSION >= 32
+/**
+   C++ Interface for Homography estimation
+
+   \param[out] H is a 3x3 array containing the estimated homography.
+   \param[out] inliers is the number of inliers that the homography was estimated to comprise,
+               in the case that htype is AF_HOMOGRAPHY_RANSAC, a higher inlier_thr value will increase the
+               estimated inliers. Note that if the number of inliers is too low, it is likely
+               that a bad homography will be returned.
+   \param[in]  x_src x coordinates of the source points.
+   \param[in]  y_src y coordinates of the source points.
+   \param[in]  x_dst x coordinates of the destination points.
+   \param[in]  y_dst y coordinates of the destination points.
+   \param[in]  htype can be AF_HOMOGRAPHY_RANSAC, for which a RANdom SAmple Consensus will be
+               used to evaluate the homography quality (e.g., number of inliers), or AF_HOMOGRAPHY_LMEDS,
+               which will use Least Median of Squares method to evaluate homography quality
+   \param[in]  inlier_thr if htype is AF_HOMOGRAPHY_RANSAC, this parameter will five the maximum L2-distance
+               for a point to be considered an inlier.
+   \param[in]  iterations maximum number of iterations when htype is AF_HOMOGRAPHY_RANSAC and backend is CPU,
+               if backend is CUDA or OpenCL, iterations is the total number of iterations, an
+               iteration is a selection of 4 random points for which the homography is estimated
+               and evaluated for number of inliers.
+   \param[in]  otype the array type for the homography output.
+
+   \ingroup cv_func_homography
+*/
+AFAPI void homography(array& H, int& inliers, const array& x_src, const array& y_src,
+                      const array& x_dst, const array& y_dst, const af_homography_type htype=AF_HOMOGRAPHY_RANSAC,
+                      const float inlier_thr=3.f, const unsigned iterations=1000, const dtype otype=f32);
+#endif
+
 }
 #endif
 
@@ -312,7 +356,8 @@ extern "C" {
 
         \ingroup cv_func_fast
     */
-    AFAPI af_err af_fast(af_features *out, const af_array in, const float thr, const unsigned arc_length, const bool non_max, const float feature_ratio, const unsigned edge);
+    AFAPI af_err af_fast(af_features *out, const af_array in, const float thr, const unsigned arc_length,
+                         const bool non_max, const float feature_ratio, const unsigned edge);
 
 #if AF_API_VERSION >= 31
     /**
@@ -341,7 +386,9 @@ extern "C" {
 
         \ingroup cv_func_harris
     */
-    AFAPI af_err af_harris(af_features *out, const af_array in, const unsigned max_corners, const float min_response, const float sigma, const unsigned block_size, const float k_thr);
+    AFAPI af_err af_harris(af_features *out, const af_array in, const unsigned max_corners,
+                           const float min_response, const float sigma,
+                           const unsigned block_size, const float k_thr);
 #endif
 
     /**
@@ -366,7 +413,9 @@ extern "C" {
 
         \ingroup cv_func_orb
     */
-    AFAPI af_err af_orb(af_features *feat, af_array *desc, const af_array in, const float fast_thr, const unsigned max_feat, const float scl_fctr, const unsigned levels, const bool blur_img);
+    AFAPI af_err af_orb(af_features *feat, af_array *desc, const af_array in,
+                        const float fast_thr, const unsigned max_feat, const float scl_fctr,
+                        const unsigned levels, const bool blur_img);
 
 #if AF_API_VERSION >= 31
     /**
@@ -400,7 +449,10 @@ extern "C" {
 
         \ingroup cv_func_sift
     */
-    AFAPI af_err af_sift(af_features *feat, af_array *desc, const af_array in, const unsigned n_layers, const float contrast_thr, const float edge_thr, const float init_sigma, const bool double_input, const float intensity_scale, const float feature_ratio);
+    AFAPI af_err af_sift(af_features *feat, af_array *desc, const af_array in,
+                         const unsigned n_layers, const float contrast_thr, const float edge_thr,
+                         const float init_sigma, const bool double_input,
+                         const float intensity_scale, const float feature_ratio);
 #endif
 
 #if AF_API_VERSION >= 32
@@ -435,7 +487,10 @@ extern "C" {
 
         \ingroup cv_func_sift
     */
-    AFAPI af_err af_gloh(af_features *feat, af_array *desc, const af_array in, const unsigned n_layers, const float contrast_thr, const float edge_thr, const float init_sigma, const bool double_input, const float intensity_scale, const float feature_ratio);
+    AFAPI af_err af_gloh(af_features *feat, af_array *desc, const af_array in,
+                         const unsigned n_layers, const float contrast_thr,
+                         const float edge_thr, const float init_sigma, const bool double_input,
+                         const float intensity_scale, const float feature_ratio);
 #endif
 
     /**
@@ -511,7 +566,8 @@ extern "C" {
 
        \ingroup cv_func_match_template
     */
-    AFAPI af_err af_match_template(af_array *out, const af_array search_img, const af_array template_img, const af_match_type m_type);
+    AFAPI af_err af_match_template(af_array *out, const af_array search_img,
+                                   const af_array template_img, const af_match_type m_type);
 
 #if AF_API_VERSION >= 31
     /**
@@ -532,7 +588,8 @@ extern "C" {
 
        \ingroup cv_func_susan
     */
-    AFAPI af_err af_susan(af_features* out, const af_array in, const unsigned radius, const float diff_thr, const float geom_thr,
+    AFAPI af_err af_susan(af_features* out, const af_array in, const unsigned radius,
+                          const float diff_thr, const float geom_thr,
                           const float feature_ratio, const unsigned edge);
 #endif
 
@@ -552,6 +609,40 @@ extern "C" {
     AFAPI af_err af_dog(af_array *out, const af_array in, const int radius1, const int radius2);
 #endif
 
+#if AF_API_VERSION >= 32
+    /**
+       C Interface wrapper for Homography estimation
+
+       \param[out] H is a 3x3 array containing the estimated homography.
+       \param[out] inliers is the number of inliers that the homography was estimated to comprise,
+                   in the case that htype is AF_HOMOGRAPHY_RANSAC, a higher inlier_thr value will increase the
+                   estimated inliers. Note that if the number of inliers is too low, it is likely
+                   that a bad homography will be returned.
+       \param[in]  x_src x coordinates of the source points.
+       \param[in]  y_src y coordinates of the source points.
+       \param[in]  x_dst x coordinates of the destination points.
+       \param[in]  y_dst y coordinates of the destination points.
+       \param[in]  htype can be AF_HOMOGRAPHY_RANSAC, for which a RANdom SAmple Consensus will be
+                   used to evaluate the homography quality (e.g., number of inliers), or AF_HOMOGRAPHY_LMEDS,
+                   which will use Least Median of Squares method to evaluate homography quality.
+       \param[in]  inlier_thr if htype is AF_HOMOGRAPHY_RANSAC, this parameter will five the maximum L2-distance
+                   for a point to be considered an inlier.
+       \param[in]  iterations maximum number of iterations when htype is AF_HOMOGRAPHY_RANSAC and backend is CPU,
+                   if backend is CUDA or OpenCL, iterations is the total number of iterations, an
+                   iteration is a selection of 4 random points for which the homography is estimated
+                   and evaluated for number of inliers.
+       \param[in]  otype the array type for the homography output.
+       \return     \ref AF_SUCCESS if the computation is is successful,
+                   otherwise an appropriate error code is returned.
+
+       \ingroup cv_func_homography
+     */
+    AFAPI af_err af_homography(af_array *H, int *inliers, const af_array x_src, const af_array y_src,
+                               const af_array x_dst, const af_array y_dst,
+                               const af_homography_type htype, const float inlier_thr,
+                               const unsigned iterations, const af_dtype otype);
+#endif
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/api/c/homography.cpp b/src/api/c/homography.cpp
@@ -0,0 +1,88 @@
+/*******************************************************
+ * Copyright (c) 2015, ArrayFire
+ * All rights reserved.
+ *
+ * This file is distributed under 3-clause BSD license.
+ * The complete license agreement can be obtained at:
+ * http://arrayfire.com/licenses/BSD-3-Clause
+ ********************************************************/
+
+#include <af/array.h>
+#include <af/defines.h>
+#include <af/vision.h>
+#include <err_common.hpp>
+#include <handle.hpp>
+#include <backend.hpp>
+#include <ArrayInfo.hpp>
+#include <homography.hpp>
+
+using af::dim4;
+using namespace detail;
+
+template<typename T>
+static inline void homography(af_array &H, int &inliers,
+                              const af_array x_src, const af_array y_src,
+                              const af_array x_dst, const af_array y_dst,
+                              const af_homography_type htype, const float inlier_thr,
+                              const unsigned iterations)
+{
+    Array<T> bestH = createEmptyArray<T>(af::dim4(3, 3));
+
+    inliers = homography<T>(bestH,
+                            getArray<float>(x_src), getArray<float>(y_src),
+                            getArray<float>(x_dst), getArray<float>(y_dst),
+                            htype, inlier_thr, iterations);
+
+    H = getHandle<T>(bestH);
+}
+
+af_err af_homography(af_array *H, int *inliers,
+                     const af_array x_src, const af_array y_src,
+                     const af_array x_dst, const af_array y_dst,
+                     const af_homography_type htype, const float inlier_thr,
+                     const unsigned iterations, const af_dtype otype)
+{
+    try {
+        ArrayInfo xsinfo = getInfo(x_src);
+        ArrayInfo ysinfo = getInfo(y_src);
+        ArrayInfo xdinfo = getInfo(x_dst);
+        ArrayInfo ydinfo = getInfo(y_dst);
+
+        af::dim4 xsdims  = xsinfo.dims();
+        af::dim4 ysdims  = ysinfo.dims();
+        af::dim4 xddims  = xdinfo.dims();
+        af::dim4 yddims  = ydinfo.dims();
+
+        af_dtype xstype = xsinfo.getType();
+        af_dtype ystype = ysinfo.getType();
+        af_dtype xdtype = xdinfo.getType();
+        af_dtype ydtype = ydinfo.getType();
+
+        if (xstype != f32) { TYPE_ERROR(1, xstype); }
+        if (ystype != f32) { TYPE_ERROR(2, ystype); }
+        if (xdtype != f32) { TYPE_ERROR(3, xdtype); }
+        if (ydtype != f32) { TYPE_ERROR(4, ydtype); }
+
+        ARG_ASSERT(1, (xsdims[0] > 0));
+        ARG_ASSERT(2, (ysdims[0] == xsdims[0]));
+        ARG_ASSERT(3, (xddims[0] > 0));
+        ARG_ASSERT(4, (yddims[0] == yddims[0]));
+
+        ARG_ASSERT(5, (inlier_thr >= 0.1f));
+        ARG_ASSERT(6, (iterations > 0));
+
+        af_array outH;
+        int outInl;
+
+        switch(otype) {
+            case f32: homography<float >(outH, outInl, x_src, y_src, x_dst, y_dst, htype, inlier_thr, iterations);  break;
+            case f64: homography<double>(outH, outInl, x_src, y_src, x_dst, y_dst, htype, inlier_thr, iterations);  break;
+            default:  TYPE_ERROR(1, otype);
+        }
+        std::swap(*H, outH);
+        std::swap(*inliers, outInl);
+    }
+    CATCHALL;
+
+    return AF_SUCCESS;
+}
diff --git a/src/api/cpp/homography.cpp b/src/api/cpp/homography.cpp
@@ -0,0 +1,32 @@
+/*******************************************************
+ * Copyright (c) 2015, ArrayFire
+ * All rights reserved.
+ *
+ * This file is distributed under 3-clause BSD license.
+ * The complete license agreement can be obtained at:
+ * http://arrayfire.com/licenses/BSD-3-Clause
+ ********************************************************/
+
+#include <af/vision.h>
+#include <af/array.h>
+#include "error.hpp"
+
+namespace af
+{
+
+void homography(array &H, int &inliers,
+                const array &x_src, const array &y_src,
+                const array &x_dst, const array &y_dst,
+                const af_homography_type htype, const float inlier_thr,
+                const unsigned iterations, const af::dtype otype)
+{
+    af_array outH;
+    AF_THROW(af_homography(&outH, &inliers,
+                           x_src.get(), y_src.get(),
+                           x_dst.get(), y_dst.get(),
+                           htype, inlier_thr, iterations, otype));
+
+    H = array(outH);
+}
+
+}
diff --git a/src/api/unified/vision.cpp b/src/api/unified/vision.cpp
@@ -66,3 +66,10 @@ af_err af_dog(af_array *out, const af_array in, const int radius1, const int rad
 {
     return CALL(out, in, radius1, radius2);
 }
+
+af_err af_homography(af_array *H, int *inliers, const af_array x_src, const af_array y_src,
+                     const af_array x_dst, const af_array y_dst, const af_homography_type htype,
+                     const float inlier_thr, const unsigned iterations, const af_dtype type)
+{
+    return CALL(H, inliers, x_src, y_src, x_dst, y_dst, htype, inlier_thr, iterations, type);
+}