Matcher ported, first page od documentation written.

modules/line_descriptor/doc/line_descriptor.rst

@@ -77,3 +77,55 @@ Apart from fields inspired to KeyPoint class, KeyLines stores information about
         /* constructor */
         KeyLine(){}
     };
+
+
+Lines extraction methodology
+----------------------------
+
+The lines extraction methodology described in the following is mainly based on [LBD]_.
+The extraction starts with a Gaussian pyramid generated from an original image, downsampled and blurred N-1 times, to obtain N layers (one for each octave), with layer 0 corresponding to input image. Then, from each layer (octave) in the pyramid, lines are extracted using LSD algorithm. 
+
+Differently from EDLine lines extractor used in original article, LSD furnishes information only about lines extremes; thus, additional information regarding slope and equation of line are computed via analytic methods. The number of pixels is obtained using `LineIterator <http://docs.opencv.org/modules/core/doc/drawing_functions.html#lineiterator>`_. Later on, all extracted lines are arranged in buckets: two lines fall in the same bucket if they represent the same line in different octave  (they have the same direction and belong to same region of original image). The set of buckets becomes the input for descriptors computation.
+
+
+Computation of binary descriptors
+---------------------------------
+
+To obtatin a binary descriptor representing a certain line detected from a certain octave of an image, we first compute a non-binary descriptor as described in [LBD]_. Given a line, we consider a rectangular region centered at it and called *line support region (LSR)*. Such region is divided into a set of bands :math:`\{B_1, B_2, ..., B_m\}`, whose length equals the one of line.
+
+If we indicate with :math:`\bf{d}_L` the direction of line, the orthogonal and clockwise direction to line :math:`\bf{d}_{\perp}` can be determined; these two directions, are used to construct a reference frame centered in the middle point of line. The gradients of pixels :math:`\bf{g'}` inside LSR can be projected to the newly determined frame, obtaining their local equivalent :math:`\bf{g'} = (\bf{g}^T \cdot \bf{d}_{\perp}, \bf{g}^T \cdot \bf{d}_L)^T \triangleq (\bf{g'}_{d_{\perp}}, \bf{g'}_{d_L})^T`.
+
+Later on, a Gaussian function is applied to all LSR's pixels along :math:`\bf{d}_\perp` direction; first, we assign a global weighting coefficient :math:`f_g(i) = (1/\sqrt{2\pi}\sigma_g)e^{-d^2_i/2\sigma^2_g}` to *i*-th row in LSR, where :math:`d_i` is the distance of *i*-th row from the center row in LSR, :math:`\sigma_g = 0.5(m \cdot w - 1)` and :math:`w` is the width of bands (the same for every band). Secondly, considering a band :math:`B_j` and its neighbor bands :math:`B_{j-1}, B_{j+1}`, we assign a local weighting  :math:`F_l(k) = (1/\sqrt{2\pi}\sigma_l)e^{-d'^2_k/2\sigma_l^2}`, where :math:`d'_k` is the distance of *k*-th row from the center row in :math:`B_j` and :math:`\sigma_l = w`. Using the global and local weights, we obtain, at the same time, the reduction of role played by gradients far from line and of boundary effect, respectively.
+
+Each band :math:`B_j` in LSR has an associated *band descriptor(BD)* which is computed considering previous and next band (top and bottom bands are ignored when computing descriptor for first and last band). Once each band has been assignen its BD, the LBD descriptor of line is simply given by
+
+.. math::
+	LBD = (BD_1^T, BD_2^T, ... , BD^T_m)^T.
+
+To compute a band descriptor :math:`B_j`, each *k*-th row in it is considered and the gradients in such row are accumulated:
+
+.. math::
+	\begin{matrix} \bf{V1}^k_j = \lambda \sum\limits_{\bf{g}'_{d_\perp}>0}\bf{g}'_{d_\perp}, &  \bf{V2}^k_j = \lambda \sum\limits_{\bf{g}'_{d_\perp}<0} -\bf{g}'_{d_\perp}, \\ \bf{V3}^k_j = \lambda \sum\limits_{\bf{g}'_{d_L}>0}\bf{g}'_{d_L}, & \bf{V4}^k_j = \lambda \sum\limits_{\bf{g}'_{d_L}<0} -\bf{g}'_{d_L}\end{matrix}.
+
+with :math:`\lambda = f_g(k)f_l(k)`.
+
+By stacking previous results, we obtain the *band description matrix (BDM)*
+
+.. math::
+	BDM_j = \left(\begin{matrix} \bf{V1}_j^1 & \bf{V1}_j^2 & \ldots & \bf{V1}_j^n \\ \bf{V2}_j^1 & \bf{V2}_j^2 & \ldots & \bf{V2}_j^n \\ \bf{V3}_j^1 & \bf{V3}_j^2 & \ldots & \bf{V3}_j^n \\ \bf{V4}_j^1 & \bf{V4}_j^2 & \ldots & \bf{V4}_j^n \end{matrix} \right) \in \mathbb{R}^{4\times n},
+
+with :math:`n` the number of rows in band :math:`B_j`:
+
+.. math::
+	n = \begin{cases} 2w, & j = 1||m; \\ 3w, & \mbox{else}. \end{cases}
+
+Each :math:`BD_j` can be obtained using the standard deviation vector :math:`S_j` and mean vector :math:`M_j` of :math:`BDM_J`. Thus, finally:
+
+.. math::
+	LBD = (M_1^T, S_1^T, M_2^T, S_2^T, \ldots, M_m^T, S_m^T)^T \in \mathbb{R}^{8m}
+
+
+References
+----------
+
+.. [LBD] Zhang, Lilian, and Reinhard Koch. *An efficient and robust line segment matching approach based on LBD descriptor and pairwise geometric consistency*, Journal of Visual Communication and Image Representation 24.7 (2013): 794-805.

modules/line_descriptor/doc/line_descriptor.rst~

@@ -10,7 +10,7 @@ Introduction
 
 One of the most challenging activities in computer vision is the extraction of useful information from a given image. Such information, usually comes in the form of points that preserve some kind of property (for instance, they are scale-invariant) and are actually representative of input image.
 
-The goal of this module is seeking a new kind of representative information inside an image and providing the functionalities for its extraction and representation. In particular, differently from previous methods for detection of relevant elements inside an image, lines are extracted in spite of points; a new class is defined ad hoc to summarize a line's properties, for reuse and plotting purposes.
+The goal of this module is seeking a new kind of representative information inside an image and providing the functionalities for its extraction and representation. In particular, differently from previous methods for detection of relevant elements inside an image, lines are extracted in place of points; a new class is defined ad hoc to summarize a line's properties, for reuse and plotting purposes.
 
 
 A class to represent a line: KeyLine
@@ -77,3 +77,49 @@ Apart from fields inspired to KeyPoint class, KeyLines stores information about
         /* constructor */
         KeyLine(){}
     };
+
+
+Lines extraction methodology
+----------------------------
+
+The lines extraction methodology described in the following is mainly based on [LBD]_.
+The extraction starts with a Gaussian pyramid generated from an original image, downsampled and blurred N-1 times, to obtain N layers (one for each octave), with layer 0 corresponding to input image. Then, from each layer (octave) in the pyramid, lines are extracted using LSD algorithm. 
+
+Differently from EDLine lines extractor used in original article, LSD furnishes information only about lines extremes; thus, additional information regarding slope and equation of line are computed via analytic methods. The number of pixels is obtained using `LineIterator <http://docs.opencv.org/modules/core/doc/drawing_functions.html#lineiterator>`_. Later on, all extracted lines are arranged in buckets: two lines fall in the same bucket if they represent the same line in different octave  (they have the same direction and belong to same region of original image). The set of buckets becomes the input for descriptors computation.
+
+
+Computation of binary descriptors
+---------------------------------
+
+To obtatin a binary descriptor representing a certain line detected from a certain octave of an image, we first compute a non-binary descriptor as described in [LBD]_. Given a line, we consider a rectangular region centered at it and called *line support region (LSR)*. Such region is divided into a set of bands :math:`\{B_1, B_2, ..., B_m\}`, whose length equals the one of line.
+
+If we indicate with :math:`\bf{d}_L` the direction of line, the orthogonal and clockwise direction to line :math:`\bf{d}_{\perp}` can be determined; these two directions, are used to construct a reference frame centered in the middle point of line. The gradients of pixels :math:`\bf{g'}` inside LSR can be projected to the newly determined frame, obtaining their local equivalent :math:`\bf{g'} = (\bf{g}^T \cdot \bf{d}_{\perp}, \bf{g}^T \cdot \bf{d}_L)^T \triangleq (\bf{g'}_{d_{\perp}}, \bf{g'}_{d_L})^T`.
+
+Later on, a Gaussian function is applied to all LSR's pixels along :math:`\bf{d}_\perp` direction; first, we assign a global weighting coefficient :math:`f_g(i) = (1/\sqrt{2\pi}\sigma_g)e^{-d^2_i/2\sigma^2_g}` to *i*-th row in LSR, where :math:`d_i` is the distance of *i*-th row from the center row in LSR, :math:`\sigma_g = 0.5(m \cdot w - 1)` and :math:`w` is the width of bands (the same for every band). Secondly, considering a band :math:`B_j` and its neighbor bands :math:`B_{j-1}, B_{j+1}`, we assign a local weighting  :math:`F_l(k) = (1/\sqrt{2\pi}\sigma_l)e^{-d'^2_k/2\sigma_l^2}`, where :math:`d'_k` is the distance of *k*-th row from the center row in :math:`B_j` and :math:`\sigma_l = w`. Using the global and local weights, we obtain, at the same time, the reduction of role played by gradients far from line and of boundary effect, respectively.
+
+Each band :math:`B_j` in LSR has an associated *band descriptor(BD)* which is computed considering previous and next band (top and bottom bands are ignored when computing descriptor for first and last band). Once each band has been assignen its BD, the LBD descriptor of line is simply given by
+
+.. math::
+	LBD = (BD_1^T, BD_2^T, ... , BD^T_m)^T.
+
+To compute a band descriptor :math:`B_j`, each *k*-th row in it is considered and the gradients in such row are accumulated:
+
+.. math::
+	\begin{matrix} \bf{V1}^k_j = \lambda \sum\limits_{\bf{g}'_{d_\perp}>0}\bf{g}'_{d_\perp}, &  \bf{V2}^k_j = \lambda \sum\limits_{\bf{g}'_{d_\perp}<0} -\bf{g}'_{d_\perp}, \\ \bf{V3}^k_j = \lambda \sum\limits_{\bf{g}'_{d_L}>0}\bf{g}'_{d_L}, & \bf{V4}^k_j = \lambda \sum\limits_{\bf{g}'_{d_L}<0} -\bf{g}'_{d_L}\end{matrix}.
+
+with :math:`\lambda = f_g(k)f_l(k)`.
+
+By stacking previous results, we obtain the *band description matrix (BDM)*
+
+.. math::
+	BDM_j = \left(\begin{matrix} \bf{V1}_j^1 & \bf{V1}_j^2 & \ldots & \bf{V1}_j^n \\ \bf{V2}_j^1 & \bf{V2}_j^2 & \ldots & \bf{V2}_j^n \\ \bf{V3}_j^1 & \bf{V3}_j^2 & \ldots & \bf{V3}_j^n \\ \bf{V4}_j^1 & \bf{V4}_j^2 & \ldots & \bf{V4}_j^n \end{matrix} \right) \in \mathbb{R}^{4\times n},
+
+with :math:`n` the number of rows in band :math:`B_j`:
+
+.. math::
+	n = \begin{cases} 2w, & j = 1||m; \\ 3w, & \mbox{else}. \end{cases}
+
+Each :math:`BD_j` can be obtained using the standard deviation vector :math:`S_j` and mean vector :math:`M_j` of :math:`BDM_J`. Thus, finally:
+
+.. math::
+	LBD = (M_1^T, S_1^T, M_2^T, S_2^T, \ldots, M_m^T, S_m^T)^T \in \mathbb{R}^{8m}

modules/line_descriptor/include/opencv2/line_descriptor/array32.hpp

@@ -0,0 +1,60 @@
+/* dynamic array of 32-bit integers
+ * arr[0]   : array size
+ * arr[1]   : array capacity
+ * arr[2..] : array content */
+
+#ifndef __OPENCV_ARRAY32_HPP
+#define __OPENCV_ARRAY32_HPP
+
+#include "types.hpp"
+
+class Array32 {
+
+ private:
+    static double ARRAY_RESIZE_FACTOR;
+    static double ARRAY_RESIZE_ADD_FACTOR;
+
+ public:
+    /* set ARRAY_RESIZE_FACTOR */
+    static void setArrayResizeFactor(double arf);
+
+    /* constructor */
+    Array32();
+
+    /* destructor */
+    ~Array32();
+
+    /* cleaning function used in destructor */
+    void cleanup();
+
+    /* push data */
+    void push(UINT32 data);
+
+    /* insert data at given index */
+    void insert(UINT32 index, UINT32 data);
+
+    /* return data */
+    UINT32* data();
+
+    /* return data size */
+    UINT32 size();
+
+    /* return capacity */
+    UINT32 capacity();
+
+    /* definition of operator = */
+    void operator= (const Array32&);
+
+    /* print data */
+    void print();
+
+    /* initializer */
+    void init(int size);
+
+    /* data */
+    UINT32 *arr;
+
+
+};
+
+#endif

modules/line_descriptor/include/opencv2/line_descriptor/bitarray.hpp

@@ -0,0 +1,66 @@
+#ifndef __OPENCV_BITARRAY_HPP
+#define __OPENCV_BITARRAY_HPP
+
+#include "types.hpp"
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+
+/* class defining a sequence of bits */
+class bitarray {
+
+ public:
+    /* pointer to bits sequence and sequence's length */
+    UINT32 *arr;
+    UINT32 length;
+	
+    /* constructor setting default values */
+    bitarray()
+    {
+        arr = NULL;
+        length = 0;
+    }
+	
+    /* constructor setting sequence's length */
+    bitarray(UINT64 _bits) {
+        init(_bits);
+    }
+	
+    /* initializer of private fields */
+    void init(UINT64 _bits)
+    {
+        length = (UINT32)ceil(_bits/32.00);
+        arr = new UINT32[length];
+        erase();
+    }
+	
+    /* destructor */
+    ~bitarray() {
+        if (arr)
+            delete[] arr;
+    }
+	
+    inline void flip(UINT64 index)
+    {
+        arr[index >> 5] ^= ((UINT32)0x01) << (index % 32);
+    }
+
+    inline void set(UINT64 index)
+    {
+        arr[index >> 5] |= ((UINT32)0x01) << (index % 32);
+    }
+	
+    inline UINT8 get(UINT64 index)
+    {
+        return (arr[index >> 5] & (((UINT32)0x01) << (index % 32))) != 0;
+    }
+	
+    /* reserve menory for an UINT32 */
+    inline void erase()
+    {
+        memset(arr, 0, sizeof(UINT32) * length);
+    }
+
+};
+
+#endif

modules/line_descriptor/include/opencv2/line_descriptor/bitops.hpp

@@ -0,0 +1,135 @@
+#ifndef __OPENCV_BITOPTS_HPP
+#define __OPENCV_BITOPTS_HPP
+
+#define popcntll __builtin_popcountll
+#define popcnt __builtin_popcount
+
+#include "precomp.hpp"
+
+/* LUT */
+const int lookup [] = {0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,1,2,
+                       2,3,2,3,3,4,2,3,3,4,3,4,4,5,1,2,2,3,
+                       2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,
+                       4,5,3,4,4,5,4,5,5,6,1,2,2,3,2,3,3,4,
+                       2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,
+                       4,5,4,5,5,6,2,3,3,4,3,4,4,5,3,4,4,5,
+                       4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,
+                       6,7,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
+                       2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,2,3,
+                       3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,
+                       4,5,5,6,4,5,5,6,5,6,6,7,2,3,3,4,3,4,
+                       4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,
+                       4,5,5,6,5,6,6,7,3,4,4,5,4,5,5,6,4,5,
+                       5,6,5,6,6,7,4,5,5,6,5,6,6,7,5,6,6,7,
+                       6,7,7,8};
+
+/*matching function */
+inline int match(UINT8*P, UINT8*Q, int codelb)
+{
+    switch(codelb)
+    {
+        case 4: // 32 bit
+            return popcnt(*(UINT32*)P ^ *(UINT32*)Q);
+            break;
+        case 8: // 64 bit
+            return popcntll(((UINT64*)P)[0] ^ ((UINT64*)Q)[0]);
+            break;
+        case 16: // 128 bit
+            return popcntll(((UINT64*)P)[0] ^ ((UINT64*)Q)[0]) \
+                + popcntll(((UINT64*)P)[1] ^ ((UINT64*)Q)[1]);
+            break;
+        case 32: // 256 bit
+            return popcntll(((UINT64*)P)[0] ^ ((UINT64*)Q)[0]) \
+                + popcntll(((UINT64*)P)[1] ^ ((UINT64*)Q)[1]) \
+                + popcntll(((UINT64*)P)[2] ^ ((UINT64*)Q)[2]) \
+                + popcntll(((UINT64*)P)[3] ^ ((UINT64*)Q)[3]);
+            break;
+        case 64: // 512 bit
+            return popcntll(((UINT64*)P)[0] ^ ((UINT64*)Q)[0]) \
+                + popcntll(((UINT64*)P)[1] ^ ((UINT64*)Q)[1]) \
+                + popcntll(((UINT64*)P)[2] ^ ((UINT64*)Q)[2]) \
+                + popcntll(((UINT64*)P)[3] ^ ((UINT64*)Q)[3]) \
+                + popcntll(((UINT64*)P)[4] ^ ((UINT64*)Q)[4]) \
+                + popcntll(((UINT64*)P)[5] ^ ((UINT64*)Q)[5]) \
+                + popcntll(((UINT64*)P)[6] ^ ((UINT64*)Q)[6]) \
+                + popcntll(((UINT64*)P)[7] ^ ((UINT64*)Q)[7]);
+            break;
+        default:
+            int output = 0;
+            for (int i=0; i<codelb; i++)
+                output+= lookup[P[i] ^ Q[i]];
+            return output;
+            break;
+    }
+
+    return -1;
+}
+
+/* splitting function (b <= 64) */
+inline void split (UINT64 *chunks, UINT8 *code, int m, int mplus, int b)
+{
+    UINT64 temp = 0x0;
+    int nbits = 0;
+    int nbyte = 0;
+    UINT64 mask = b==64 ? 0xFFFFFFFFFFFFFFFFLLU : ((UINT64_1 << b) - UINT64_1);
+
+    for (int i=0; i<m; i++)
+    {
+        while (nbits < b)
+        {
+            temp |= ((UINT64)code[nbyte++] << nbits);
+            nbits += 8;
+        }
+
+        chunks[i] = temp & mask;
+        temp = b==64 ? 0x0 : temp >> b;
+        nbits -= b;
+
+        if (i == mplus-1)
+        {
+            b--; /* b <= 63 */
+            mask = ((UINT64_1 << b) - UINT64_1);
+        }
+    }
+}
+
+/* generates the next binary code (in alphabetical order) with the
+   same number of ones as the input x. Taken from
+   http://www.geeksforgeeks.org/archives/10375 */
+inline UINT64 next_set_of_n_elements(UINT64 x)
+{
+    UINT64 smallest, ripple, new_smallest;
+
+    smallest     = x & -x;
+    ripple       = x + smallest;
+    new_smallest = x ^ ripple;
+    new_smallest = new_smallest / smallest;
+    new_smallest >>= 2;
+    return ripple | new_smallest;
+}
+
+/* print code */
+inline void print_code(UINT64 tmp, int b)
+{
+    for (int j=(b-1); j>=0; j--)
+    {
+        printf("%llu", (long long int) tmp/(1 << j));
+        tmp = tmp - (tmp/(1 << j)) * (1 << j);
+    }
+
+    printf("\n");
+}
+
+inline UINT64 choose(int n, int r)
+{
+    UINT64 nchooser = 1;
+    for (int k=0; k < r; k++)
+    {
+        nchooser *= n-k;
+        nchooser /= k+1;
+    }
+
+    return nchooser;
+}
+
+#endif

modules/line_descriptor/include/opencv2/line_descriptor/bucket_group.hpp

@@ -0,0 +1,29 @@
+#ifndef __OPENCV_BUCKET_GROUP_HPP
+#define __OPENCV_BUCKET_GROUP_HPP
+
+#include "types.hpp"
+#include "array32.hpp"
+#include "bitarray.hpp"
+
+class BucketGroup {
+
+ public:
+    /* constructor */
+    BucketGroup();
+
+    /* destructor */
+    ~BucketGroup();
+
+    /* insert data into the bucket */
+    void insert(int subindex, UINT32 data);
+
+    /* perform a query to the bucket */
+    UINT32* query(int subindex, int *size);
+
+    /* data fields */
+    UINT32 empty;
+    Array32 *group;
+
+};
+
+#endif

modules/line_descriptor/include/opencv2/line_descriptor/descriptor.hpp

@@ -47,9 +47,13 @@
 #define __OPENCV_DESCRIPTOR_HPP__
 
 #include "LineStructure.hpp"
-#include "opencv2/core.hpp"
+#include "array32.hpp"
-#include <opencv2/features2d.hpp>
+#include "bitarray.hpp"
-
+#include "bitops.hpp"
+#include "bucket_group.hpp"
+#include "mihasher.hpp"
+#include "sparse_hashtable.hpp"
+#include "types.hpp"
 
 
 namespace cv
@@ -107,18 +111,9 @@ namespace cv
         struct CV_EXPORTS_W_SIMPLE Params{
             CV_WRAP Params();
 
-            /* global threshold for line descriptor distance, default is 0.35 */
-            CV_PROP_RW  float LowestThreshold;
-
-            /* the NNDR threshold for line descriptor distance, default is 0.6 */
-            CV_PROP_RW float NNDRThreshold;
-
             /* the number of image octaves (default = 5) */
             CV_PROP_RW int  numOfOctave_;
 
-            /* the number of bands used to compute line descriptor (default: 9) */
-            CV_PROP_RW int  numOfBand_;
-
             /* the width of band; (default: 7) */
             CV_PROP_RW int  widthOfBand_;
 
@@ -144,6 +139,14 @@ namespace cv
         /* destructor */
         ~BinaryDescriptor();
 
+        /* setters and getters */
+        int getNumOfOctaves();
+        void setNumOfOctaves(int octaves);
+        int getWidthOfBand();
+        void setWidthOfBand(int width);
+        int getReductionRatio();
+        void setReductionRatio(int rRatio);
+
         /* read parameters from a FileNode object and store them (class function ) */
         virtual void read( const cv::FileNode& fn );
 
@@ -248,6 +251,76 @@ namespace cv
 
     };
 
+    class CV_EXPORTS_W BinaryDescriptorMatcher: public Algorithm
+    {
+
+    public:
+        /* for every input descriptor,
+           find the best matching one (for a pair of images) */
+        void match( const Mat& queryDescriptors,
+                    const Mat& trainDescriptors,
+                    std::vector<DMatch>& matches,
+                    const Mat& mask=Mat() ) const;
+
+        /* for every input descriptor,
+           find the best matching one (from one image to a set) */
+        void match( const Mat& queryDescriptors,
+                    std::vector<DMatch>& matches,
+                    const std::vector<Mat>& masks=std::vector<Mat>() );
+
+        /* for every input descriptor,
+           find the best k matching descriptors (for a pair of images) */
+        void knnMatch( const Mat& queryDescriptors,
+                       const Mat& trainDescriptors,
+                       std::vector<std::vector<DMatch> >& matches,
+                       int k,
+                       const Mat& mask=Mat(),
+                       bool compactResult=false ) const;
+
+        /* for every input descriptor,
+           find the best k matching descriptors (from one image to a set) */
+        void knnMatch( const Mat& queryDescriptors,
+                       std::vector<std::vector<DMatch> >& matches,
+                       int k,
+                       const std::vector<Mat>& masks=std::vector<Mat>(),
+                       bool compactResult=false );
+
+        /* for every input desciptor, find all the ones falling in a
+           certaing atching radius (for a pair of images) */
+        void radiusMatch( const Mat& queryDescriptors,
+                          const Mat& trainDescriptors,
+                          std::vector<std::vector<DMatch> >& matches,
+                          float maxDistance,
+                          const Mat& mask=Mat(),
+                          bool compactResult=false ) const;
+
+        /* for every input desciptor, find all the ones falling in a
+           certaing atching radius (from one image to a set) */
+        void radiusMatch( const Mat& queryDescriptors,
+                          std::vector<std::vector<DMatch> >& matches,
+                          float maxDistance,
+                          const std::vector<Mat>& masks=std::vector<Mat>(),
+                          bool compactResult=false );
+
+        /* constructor with smart pointer */
+        static Ptr<BinaryDescriptorMatcher> createBinaryDescriptorMatcher();
+
+
+        /* write/read data to/from file */
+        virtual void read( const FileNode& );
+        virtual void write( FileStorage& ) const;
+
+        /* constructor */
+        BinaryDescriptorMatcher(){};
+
+        /* desctructor */
+        ~BinaryDescriptorMatcher(){};
+
+    private:
+        /* vector to store new desciptors */
+        std::vector<Mat> descriptorsVector;
+    };
+
 }
 
 #endif

modules/line_descriptor/include/opencv2/line_descriptor/mihasher.hpp

@@ -0,0 +1,88 @@
+#ifndef __OPENCV_MIHASHER_HPP
+#define __OPENCV_MIHASHER_HPP
+
+#include "types.hpp"
+#include "bitops.hpp"
+
+#include "sparse_hashtable.hpp"
+#include "bitarray.hpp"
+
+#include <opencv2/core.hpp>
+#include <opencv2/imgproc.hpp>
+
+//#include <iostream>
+//#include <math.h>
+
+
+class Mihasher {
+ private:
+
+    /* Bits per code */
+    int B;
+
+    /* B/8 */
+    int B_over_8;
+
+    /* Bits per chunk (must be less than 64) */
+    int b;
+
+    /* Number of chunks */
+    int m;
+
+    /* Number of chunks with b bits (have 1 bit more than others) */
+    int mplus;
+
+    /* Maximum hamming search radius (we use B/2 by default) */
+    int D;
+
+    /* Maximum hamming search radius per substring */
+    int d;
+
+    /* Maximum results to return */
+    int K;
+
+    /* Number of codes */
+    UINT64 N;
+    
+    /* Table of original full-length codes */
+    cv::Mat codes;
+
+    /* Counter for eliminating duplicate results (it is not thread safe) */
+    bitarray *counter;
+    
+    /* Array of m hashtables */
+    SparseHashtable *H;
+        
+    /* Volume of a b-bit Hamming ball with radius s (for s = 0 to d) */
+    UINT32 *xornum;
+
+    /* Used within generation of binary codes at a certain Hamming distance */
+    int power[100];
+
+ public:
+    
+    /* constructor */
+    Mihasher();
+
+    /* desctructor */
+    ~Mihasher();
+
+    /* constructor 2 */
+    Mihasher(int B, int m);
+
+    /* K setter */
+    void setK(int K);
+
+    /* populate tables */
+    void populate(cv::Mat & codes, UINT32 N, int dim1codes);
+
+    /* execute a batch query */
+    void batchquery (UINT32 * results, UINT32 *numres/*, qstat *stats*/,const cv::Mat & q, UINT32 numq, int dim1queries);
+    
+ private:
+
+    /* execute a single query */
+    void query(UINT32 * results, UINT32* numres/*, qstat *stats*/, UINT8 *q, UINT64 * chunks, UINT32 * res, int query_i);
+};
+
+#endif

modules/line_descriptor/include/opencv2/line_descriptor/sparse_hashtable.hpp

@@ -0,0 +1,43 @@
+#ifndef __OPENCV_SPARSE_HASHTABLE_HPP
+#define __OPENCV_SPARSE_HASHTABLE_HPP
+
+#include "types.hpp"
+#include "bucket_group.hpp"
+
+class SparseHashtable
+{
+
+ private:
+
+    /* Maximum bits per key before folding the table */
+    static const int MAX_B;
+
+    /* Bins (each bin is an Array object for duplicates of the same key) */
+    BucketGroup *table;
+
+ public:
+
+    /* constructor */
+    SparseHashtable();
+
+    /* destructor */
+    ~SparseHashtable();
+
+    /* initializer */
+    int init(int _b);
+	
+    /* insert data */
+    void insert(UINT64 index, UINT32 data);
+
+    /* query data */
+    UINT32* query(UINT64 index, int* size);
+
+    /* Bits per index */
+    int b;
+
+    /*  Number of bins */
+    UINT64 size;
+
+};
+
+#endif

modules/line_descriptor/include/opencv2/line_descriptor/types.hpp

@@ -0,0 +1,16 @@
+#include <inttypes.h>
+
+#ifndef __OPENCV_TYPES_HPP
+#define __OPENCV_TYPES_HPP
+
+/* define data types */
+typedef uint64_t UINT64;
+typedef uint32_t UINT32;
+typedef uint16_t UINT16;
+typedef uint8_t UINT8;
+
+/* define constants */
+#define UINT64_1 ((UINT64)0x01)
+#define UINT32_1 ((UINT32)0x01)
+
+#endif

modules/line_descriptor/samples/compute_descriptors.cpp

@@ -10,7 +10,6 @@
 #include <iostream>
 
 using namespace cv;
-using namespace std;
 
 static const char* keys =
 {
@@ -19,10 +18,11 @@ static const char* keys =
 
 static void help()
 {
-  cout << "\nThis example shows the functionalities of lines extraction " <<
+  std::cout << "\nThis example shows the functionalities of lines extraction " <<
           "and descriptors computation furnished by BinaryDescriptor class\n" <<
           "Please, run this sample using a command in the form\n" <<
-          "./example_line_descriptor_compute_descriptors <path_to_input_image>" << endl;
+          "./example_line_descriptor_compute_descriptors <path_to_input_image>"
+          << std::endl;
 }
 
 inline void writeMat(cv::Mat m, std::string name, int n)

modules/line_descriptor/samples/lines_extraction.cpp

@@ -57,6 +57,7 @@ int main( int argc, char** argv )
 
     /* extract lines */
     bd->detect(imageMat, lines, mask);
+    std::cout << lines.size() << std::endl;
 
     /* draw lines extracted from octave 0 */
     cv::Mat output = imageMat.clone();

modules/line_descriptor/samples/matching.cpp

@@ -0,0 +1,145 @@
+#include <opencv2/line_descriptor.hpp>
+
+#include "opencv2/core/utility.hpp"
+#include "opencv2/core/private.hpp"
+#include <opencv2/imgproc.hpp>
+#include <opencv2/features2d.hpp>
+#include <opencv2/highgui.hpp>
+
+#include <iostream>
+
+using namespace cv;
+
+static const char* keys =
+{
+    "{@image_path1 | | Image path 1 }"
+    "{@image_path2 | | Image path 2 }"
+};
+
+static void help()
+{
+  std::cout << "\nThis example shows the functionalities of lines extraction " <<
+          "and descriptors computation furnished by BinaryDescriptor class\n" <<
+          "Please, run this sample using a command in the form\n" <<
+          "./example_line_descriptor_compute_descriptors <path_to_input_image 1>"
+          << "<path_to_input_image 2>" << std::endl;
+
+}
+
+inline void writeMat(cv::Mat m, std::string name, int n)
+{
+    std::stringstream ss;
+    std::string s;
+    ss << n;
+    ss >> s;
+    std::string fileNameConf = name + s;
+    cv::FileStorage fsConf(fileNameConf, cv::FileStorage::WRITE);
+    fsConf << "m" << m;
+
+    fsConf.release();
+}
+
+int main( int argc, char** argv )
+{
+    /* get parameters from comand line */
+    CommandLineParser parser( argc, argv, keys );
+    String image_path1 = parser.get<String>( 0 );
+    String image_path2 = parser.get<String>( 1 );
+
+    if(image_path1.empty() || image_path2.empty())
+    {
+        help();
+        return -1;
+    }
+
+
+    /* load image */
+    cv::Mat imageMat1 = imread(image_path1, 0);
+    cv::Mat imageMat2 = imread(image_path2, 0);
+    if(imageMat1.data == NULL || imageMat2.data == NULL)
+    {
+        std::cout << "Error, images could not be loaded. Please, check their path"
+                  << std::endl;
+    }
+
+    /* create binary masks */
+    cv::Mat mask1 = Mat::ones(imageMat1.size(), CV_8UC1);
+    cv::Mat mask2 = Mat::ones(imageMat2.size(), CV_8UC1);
+
+    /* create a pointer to a BinaryDescriptor object with default parameters */
+    Ptr<BinaryDescriptor> bd = BinaryDescriptor::createBinaryDescriptor();
+
+    /* compute lines */
+    std::vector<KeyLine> keylines1, keylines2;
+    bd->detect(imageMat1, keylines1, mask1);
+    bd->detect(imageMat2, keylines2, mask2);
+
+    std::cout << "lines " << keylines1.size() << " " << keylines2.size()
+              << std::endl;
+
+    /* compute descriptors */
+    cv::Mat descr1, descr2;
+    bd->compute(imageMat1, keylines1, descr1);
+    bd->compute(imageMat2, keylines2, descr2);
+
+    /* create a BinaryDescriptorMatcher object */
+    Ptr<BinaryDescriptorMatcher> bdm = BinaryDescriptorMatcher::createBinaryDescriptorMatcher();
+
+    /* require match */
+    std::vector<DMatch> matches;
+    bdm->match(descr1, descr2, matches);
+    for(int x = 0; x<matches.size(); x++)
+        std::cout << matches[x].queryIdx << " " << matches[x].trainIdx << std::endl;
+
+    /* result checkout */
+    cv::Mat result(descr1.size(), CV_8UC1);
+    std::cout << "size " << descr1.rows << " " << descr1.cols
+              << " " << descr2.rows << " " << descr2.cols << std::endl;
+//    for(size_t i = 0; i<matches.size(); i++){
+//        uchar* pointer = result.ptr(i);
+//        uchar* trainPointer = descr2.ptr(matches[i].trainIdx);
+//        *pointer = *trainPointer;
+//        pointer++;
+//    }
+
+
+    /* write matrices */
+    writeMat(descr1, "descr1", 0);
+    writeMat(result, "result", 0);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

modules/line_descriptor/src/BinaryDescriptor.cpp

@@ -40,6 +40,7 @@
  //M*/
 
 #define _USE_MATH_DEFINES
+#define NUM_OF_BANDS 9
 
 #include "precomp.hpp"
 
@@ -54,22 +55,46 @@ static const int combinations [32][2] = {{0,1},{0,2},{0,3},{0,4},{0,5},{0,6},{1,
 /* return default parameters */
 BinaryDescriptor::Params::Params()
 {
-    LowestThreshold = 0.35;
-    NNDRThreshold = 0.6;
     numOfOctave_ = 1;
-    numOfBand_ = 9;
     widthOfBand_ = 7;
     reductionRatio = 2;
+}
 
+/* setters and getters */
+int BinaryDescriptor::getNumOfOctaves()
+{
+    return params.numOfOctave_;
+}
+
+void BinaryDescriptor::setNumOfOctaves(int octaves)
+{
+    params.numOfOctave_ = octaves;
+}
+
+int BinaryDescriptor::getWidthOfBand()
+{
+    return params.widthOfBand_;
+}
+
+void BinaryDescriptor::setWidthOfBand(int width)
+{
+    params.widthOfBand_ = width;
+}
+
+int BinaryDescriptor::getReductionRatio()
+{
+    return params.reductionRatio;
+}
+
+void BinaryDescriptor::setReductionRatio(int rRatio)
+{
+    params.reductionRatio = rRatio;
 }
 
 /* read parameters from a FileNode object and store them (struct function) */
 void BinaryDescriptor::Params::read(const cv::FileNode& fn )
 {
-    LowestThreshold = fn["LowestThreshold"];
-    NNDRThreshold = fn["NNDRThreshold"];
     numOfOctave_ = fn["numOfOctave_"];
-    numOfBand_ = fn["numOfBand_"];
     widthOfBand_ = fn["widthOfBand_"];
     reductionRatio = fn["reductionRatio"];
 }
@@ -77,10 +102,8 @@ void BinaryDescriptor::Params::read(const cv::FileNode& fn )
 /* store parameters to a FileStorage object (struct function) */
 void BinaryDescriptor::Params::write(cv::FileStorage& fs) const
 {
-    fs << "LowestThreshold" << LowestThreshold;
-    fs << "NNDRThreshold" << NNDRThreshold;
     fs << "numOfOctave_" << numOfOctave_;
-    fs << "numOfBand_" <<  numOfBand_;
+    fs << "numOfBand_" <<  NUM_OF_BANDS;
     fs << "widthOfBand_" << widthOfBand_;
     fs << "reductionRatio" << reductionRatio;
 }
@@ -118,15 +141,15 @@ BinaryDescriptor::BinaryDescriptor(const BinaryDescriptor::Params &parameters) :
     }
 
     /* prepare a vector for global weights F_g*/
-    gaussCoefG_.resize(params.numOfBand_*params.widthOfBand_);
+    gaussCoefG_.resize(NUM_OF_BANDS*params.widthOfBand_);
 
     /* compute center of LSR */
-    u = (params.numOfBand_*params.widthOfBand_-1)/2;
+    u = (NUM_OF_BANDS*params.widthOfBand_-1)/2;
 
     /* compute exponential part of F_g */
     sigma = u;
     invsigma2 = -1/(2*sigma*sigma);
-    for(int i=0; i<params.numOfBand_*params.widthOfBand_; i++)
+    for(int i=0; i<NUM_OF_BANDS*params.widthOfBand_; i++)
     {
         dis = i-u;
         gaussCoefG_[i] = exp(dis*dis*invsigma2);
@@ -886,8 +909,8 @@ int BinaryDescriptor::computeLBD(ScaleLines &keyLines)
     short numOfFinalLine = keyLines.size();
     float *dL = new float[2];//line direction cos(dir), sin(dir)
     float *dO = new float[2];//the clockwise orthogonal vector of line direction.
-    short heightOfLSP = params.widthOfBand_*params.numOfBand_;//the height of line support region;
+    short heightOfLSP = params.widthOfBand_*NUM_OF_BANDS;//the height of line support region;
-    short descriptor_size = params.numOfBand_ * 8;//each band, we compute the m( pgdL, ngdL,  pgdO, ngdO) and std( pgdL, ngdL,  pgdO, ngdO);
+    short descriptor_size = NUM_OF_BANDS * 8;//each band, we compute the m( pgdL, ngdL,  pgdO, ngdO) and std( pgdL, ngdL,  pgdO, ngdO);
     float pgdLRowSum;//the summation of {g_dL |g_dL>0 } for each row of the region;
     float ngdLRowSum;//the summation of {g_dL |g_dL<0 } for each row of the region;
     float pgdL2RowSum;//the summation of {g_dL^2 |g_dL>0 } for each row of the region;
@@ -897,16 +920,16 @@ int BinaryDescriptor::computeLBD(ScaleLines &keyLines)
     float pgdO2RowSum;//the summation of {g_dO^2 |g_dO>0 } for each row of the region;
     float ngdO2RowSum;//the summation of {g_dO^2 |g_dO<0 } for each row of the region;
 
-    float *pgdLBandSum  = new float[params.numOfBand_];//the summation of {g_dL |g_dL>0 } for each band of the region;
+    float *pgdLBandSum  = new float[NUM_OF_BANDS];//the summation of {g_dL |g_dL>0 } for each band of the region;
-    float *ngdLBandSum  = new float[params.numOfBand_];//the summation of {g_dL |g_dL<0 } for each band of the region;
+    float *ngdLBandSum  = new float[NUM_OF_BANDS];//the summation of {g_dL |g_dL<0 } for each band of the region;
-    float *pgdL2BandSum = new float[params.numOfBand_];//the summation of {g_dL^2 |g_dL>0 } for each band of the region;
+    float *pgdL2BandSum = new float[NUM_OF_BANDS];//the summation of {g_dL^2 |g_dL>0 } for each band of the region;
-    float *ngdL2BandSum = new float[params.numOfBand_];//the summation of {g_dL^2 |g_dL<0 } for each band of the region;
+    float *ngdL2BandSum = new float[NUM_OF_BANDS];//the summation of {g_dL^2 |g_dL<0 } for each band of the region;
-    float *pgdOBandSum  = new float[params.numOfBand_];//the summation of {g_dO |g_dO>0 } for each band of the region;
+    float *pgdOBandSum  = new float[NUM_OF_BANDS];//the summation of {g_dO |g_dO>0 } for each band of the region;
-    float *ngdOBandSum  = new float[params.numOfBand_];//the summation of {g_dO |g_dO<0 } for each band of the region;
+    float *ngdOBandSum  = new float[NUM_OF_BANDS];//the summation of {g_dO |g_dO<0 } for each band of the region;
-    float *pgdO2BandSum = new float[params.numOfBand_];//the summation of {g_dO^2 |g_dO>0 } for each band of the region;
+    float *pgdO2BandSum = new float[NUM_OF_BANDS];//the summation of {g_dO^2 |g_dO>0 } for each band of the region;
-    float *ngdO2BandSum = new float[params.numOfBand_];//the summation of {g_dO^2 |g_dO<0 } for each band of the region;
+    float *ngdO2BandSum = new float[NUM_OF_BANDS];//the summation of {g_dO^2 |g_dO<0 } for each band of the region;
 
-    short numOfBitsBand = params.numOfBand_*sizeof(float);
+    short numOfBitsBand = NUM_OF_BANDS*sizeof(float);
     short lengthOfLSP; //the length of line support region, varies with lines
     short halfHeight = (heightOfLSP-1)/2;
     short halfWidth;
@@ -1075,7 +1098,7 @@ int BinaryDescriptor::computeLBD(ScaleLines &keyLines)
                 }
 
                 bandID = bandID+2;
-                if(bandID<params.numOfBand_){/*the band below the current band */
+                if(bandID<NUM_OF_BANDS){/*the band below the current band */
                     coefInGaussion = gaussCoefL_[hID%params.widthOfBand_];
                     pgdLBandSum[bandID] +=  coefInGaussion * pgdLRowSum;
                     ngdLBandSum[bandID] +=  coefInGaussion * ngdLRowSum;
@@ -1101,8 +1124,8 @@ int BinaryDescriptor::computeLBD(ScaleLines &keyLines)
             float invN2 = 1.0/(params.widthOfBand_ * 2.0);
             float invN3 = 1.0/(params.widthOfBand_ * 3.0);
             float invN, temp;
-            for(bandID = 0; bandID<params.numOfBand_; bandID++){
+            for(bandID = 0; bandID<NUM_OF_BANDS; bandID++){
-                if(bandID==0||bandID==params.numOfBand_-1){	invN = invN2;
+                if(bandID==0||bandID==NUM_OF_BANDS-1){	invN = invN2;
                 }else{ invN = invN3;}
                 desID = bandID * 8;
                 temp = pgdLBandSum[bandID] * invN;
@@ -1127,7 +1150,7 @@ int BinaryDescriptor::computeLBD(ScaleLines &keyLines)
             desVec = pSingleLine->descriptor.data();
 
             int base = 0;
-            for(short i=0; i<params.numOfBand_*8; ++base, i=base*8){
+            for(short i=0; i<NUM_OF_BANDS*8; ++base, i=base*8){
                 tempM += *(desVec+i) * *(desVec+i);//desVec[8*i+0] * desVec[8*i+0];
                 tempM += *(desVec+i+1) * *(desVec+i+1);//desVec[8*i+1] * desVec[8*i+1];
                 tempM += *(desVec+i+2) * *(desVec+i+2);//desVec[8*i+2] * desVec[8*i+2];
@@ -1142,7 +1165,7 @@ int BinaryDescriptor::computeLBD(ScaleLines &keyLines)
             tempS = 1/sqrt(tempS);
             desVec = pSingleLine->descriptor.data();
             base = 0;
-            for(short i=0; i<params.numOfBand_*8; ++base, i=base*8){
+            for(short i=0; i<NUM_OF_BANDS*8; ++base, i=base*8){
                 *(desVec+i) = *(desVec+i) * tempM;//desVec[8*i] =  desVec[8*i] * tempM;
                 *(desVec+1+i) = *(desVec+1+i) * tempM;//desVec[8*i+1] =  desVec[8*i+1] * tempM;
                 *(desVec+2+i) = *(desVec+2+i) * tempM;//desVec[8*i+2] =  desVec[8*i+2] * tempM;

modules/line_descriptor/src/BinaryDescriptorMatcher.cpp

@@ -0,0 +1,57 @@
+#include "precomp.hpp"
+
+using namespace cv;
+
+/* constructor with smart pointer */
+Ptr<BinaryDescriptorMatcher> BinaryDescriptorMatcher::createBinaryDescriptorMatcher()
+{
+    return Ptr<BinaryDescriptorMatcher>(new BinaryDescriptorMatcher());
+}
+
+void BinaryDescriptorMatcher::read( const FileNode& ){}
+void BinaryDescriptorMatcher::write( FileStorage& ) const{}
+
+/* for every input descriptor, find the best matching one (for a pair of images) */
+void BinaryDescriptorMatcher::match( const Mat& queryDescriptors,
+                                     const Mat& trainDescriptors,
+                                     std::vector<DMatch>& matches,
+                                     const Mat& mask ) const
+{
+    /* create a new mihasher object */
+    Mihasher *mh = new Mihasher(256, 32);
+
+    /* populate mihasher */
+    cv::Mat copy = trainDescriptors.clone();
+    mh->populate(copy, copy.rows, copy.cols);
+    mh->setK(1);
+
+    /* prepare structures for query */
+    UINT32 *results = new UINT32[queryDescriptors.rows];
+    UINT32 * numres = new UINT32[(256+1)*(queryDescriptors.rows)];
+
+    /* execute query */
+    mh->batchquery(results,
+                   numres,
+                   queryDescriptors,
+                   queryDescriptors.rows,
+                   queryDescriptors.cols);
+
+    /* compose matches */
+    for(size_t counter = 0; counter<queryDescriptors.rows; counter++)
+    {
+        /* create a DMatch object if required by mask of if there is
+           no mask at all */
+        if( mask.empty() || (!mask.empty() && mask.at<int>(counter)!=0))
+        {
+            DMatch dm;
+            dm.queryIdx = counter;
+            dm.trainIdx = results[counter];
+            dm.imgIdx = 0;
+            dm.distance = numres[counter];
+
+            matches.push_back(dm);
+        }
+    }
+
+
+}

modules/line_descriptor/src/array32.cpp

@@ -0,0 +1,155 @@
+/* dynamic array of 32-bit integers
+ * arr[0]   : array size
+ * arr[1]   : array capacity
+ * arr[2..] : array content */
+
+#include "precomp.hpp"
+
+
+/* no need for the static keyword in the definition */
+double Array32::ARRAY_RESIZE_FACTOR = 1.1;    // minimum is 1.0
+double Array32::ARRAY_RESIZE_ADD_FACTOR = 4;  // minimum is 1
+
+/* set ARRAY_RESIZE_FACTOR */
+void Array32::setArrayResizeFactor(double arf)
+{
+    ARRAY_RESIZE_FACTOR = arf;
+}
+
+/* constructor */
+Array32::Array32 ()
+{
+    arr = NULL;
+}
+
+/* definition of operator =
+  Array32& Array32::operator = (const Array32 &rhs) */
+void Array32::operator = (const Array32 &rhs)
+{
+    if (&rhs != this)
+    this->arr = rhs.arr;
+}
+
+/* destructor */
+Array32::~Array32 ()
+{
+    cleanup();
+}
+
+/* cleaning function used in destructor */
+void Array32::cleanup ()
+{
+    free(arr);
+}
+
+/* push data */
+void Array32::push(UINT32 Data)
+{
+    if (arr)
+    {
+    if (arr[0] == arr[1])
+    {
+        arr[1] = std::max(ceil(arr[1]*ARRAY_RESIZE_FACTOR),
+                          arr[1]+ARRAY_RESIZE_ADD_FACTOR);
+        UINT32* new_Data = static_cast<UINT32*>
+                           (realloc (arr, sizeof(UINT32)*(2 + arr[1])));
+        if (new_Data == NULL)
+        {
+            /* could not realloc, but orig still valid */
+            std::cout << "ALERT!!!! Not enough memory, operation aborted!"
+                      << std::endl;
+            exit(0);
+        }
+	     else
+	     {
+             arr = new_Data;
+	     }
+	    
+	}
+
+    arr[2 + arr[0]] = Data;
+	arr[0]++;
+
+    }
+
+    else
+    {
+        arr = (UINT32*) malloc ((2+ARRAY_RESIZE_ADD_FACTOR)*sizeof(UINT32));
+        arr[0] = 1;
+        arr[1] = 1;
+        arr[2] = Data;
+    }
+}
+
+/* insert data at given index */
+void Array32::insert(UINT32 index, UINT32 Data) {
+    if (arr) {
+	if (arr[0] == arr[1]) {
+	    arr[1] = ceil(arr[0]*1.1);
+           UINT32* new_data = static_cast<UINT32*>
+                   (realloc (arr, sizeof(UINT32)*(2 + arr[1])));
+	         if (new_data == NULL)
+	         {
+	           // could not realloc, but orig still valid
+                 std::cout << "ALERT!!!! Not enough memory, operation aborted!"
+                           << std::endl;
+	             exit(0);
+	         }
+	         else
+	         {
+	             arr = new_data;
+	         }
+	}
+
+	memmove(arr+(2+index)+1, arr+(2+index), (arr[0]-index)*sizeof(*arr));
+
+    arr[2+index] = Data;
+	arr[0]++;
+    }
+
+    else
+    {
+        arr = (UINT32*) malloc (3*sizeof(UINT32));
+        arr[0] = 1;
+        arr[1] = 1;
+        arr[2] = Data;
+    }
+}
+
+
+/* return data */
+UINT32* Array32::data()
+{
+    return arr? arr + 2 : NULL;
+}
+
+/* return data size */
+UINT32 Array32::size ()
+{
+    return arr ? arr[0] : 0;
+}
+
+/* return capacity */
+UINT32 Array32::capacity ()
+{
+    return arr ? arr[1] : 0;
+}
+
+/* print data */
+void Array32::print() {
+    for (int i=0; i<size(); i++)
+        printf("%d, ", arr[i+2]);
+
+    printf("\n");
+}
+
+/* initializer */
+void Array32::init(int size)
+{
+    if (arr == NULL)
+    {
+        arr = (UINT32*) malloc ((2+size)*sizeof(UINT32));
+        arr[0] = 0;
+        arr[1] = size;
+    }
+}

modules/line_descriptor/src/bucket_group.cpp

@@ -0,0 +1,58 @@
+#include "precomp.hpp"
+
+/* constructor */
+BucketGroup::BucketGroup()
+{
+    empty = 0;
+    group = NULL;
+}
+
+/* destructor */
+BucketGroup::~BucketGroup()
+{
+    if (group != NULL)
+        delete group;
+}
+
+/* insert data into the bucket */
+void BucketGroup::insert(int subindex, UINT32 data)
+{
+    if (group == NULL)
+    {
+        group = new Array32();
+        group->push(0);
+    }
+
+    UINT32 lowerbits = ((UINT32)1 << subindex) - 1;
+    int end = popcnt(empty & lowerbits);
+
+    if (!(empty & ((UINT32)1 << subindex)))
+    {
+        group->insert(end, group->arr[end+2]);
+        empty |= (UINT32)1 << subindex;
+    }
+
+    int totones = popcnt(empty);
+    group->insert(totones+1+group->arr[2+end+1], data);
+    for (int i=end+1; i<totones+1; i++)
+        group->arr[2+i]++;
+}
+
+/* perform a query to the bucket */
+UINT32* BucketGroup::query(int subindex, int *size)
+{
+    if (empty & ((UINT32)1 << subindex))
+    {
+        UINT32 lowerbits = ((UINT32)1 << subindex) - 1;
+        int end = popcnt(empty & lowerbits);
+        int totones = popcnt(empty);
+        *size = group->arr[2+end+1]-group->arr[2+end];
+        return group->arr + 2 + totones+1 + group->arr[2+end];
+    }
+
+    else
+    {
+        *size = 0;
+        return NULL;
+    }
+}

modules/line_descriptor/src/mihasher.cpp

@@ -0,0 +1,227 @@
+#include "precomp.hpp"
+
+
+/* execute a batch query */
+void Mihasher::batchquery(UINT32 * results, UINT32 *numres, const cv::Mat & queries, UINT32 numq, int dim1queries)
+{
+    /* create and initialize a bitarray */
+    counter = new bitarray;
+    counter->init(N);
+
+    UINT32 *res  = new UINT32[K*(D+1)];
+    UINT64 *chunks = new UINT64[m];
+    UINT32 * presults = results;
+    UINT32 *pnumres = numres;
+
+    /* make a copy of input queries */
+    cv::Mat queries_clone = queries.clone();
+
+    /* set a pointer to first query (row) */
+    UINT8 *pq = queries_clone.ptr();
+
+    /* loop over number of descriptors */
+    for (size_t i=0; i<numq; i++) {
+        /* for every descriptor, query database */
+        query(presults, pnumres, pq, chunks, res, i);
+
+        /* move pointer to write next K indeces */
+        presults += K;
+        pnumres += B+1;
+
+        /* move forward pointer to current row in descriptors matrix */
+        pq += dim1queries;
+
+    }
+
+    delete [] res;
+    delete [] chunks;
+
+    delete counter;
+}
+
+/* execute a single query */
+void Mihasher::query(UINT32*  results, UINT32* numres,
+                     UINT8 * Query, UINT64 *chunks, UINT32 *res, int query_i)
+{ 
+    /* if K == 0 that means we want everything to be processed.
+       So maxres = N in that case. Otherwise K limits the results processed */
+    UINT32 maxres = K ? K : N;
+
+    /* number of results so far obtained (up to a distance of s per chunk) */
+    UINT32 n = 0;
+
+    /* number of candidates tested with full codes (not counting duplicates) */
+    UINT32 nc = 0;
+
+    /* counting everything retrieved (duplicates are counted multiple times)
+       number of lookups (and xors) */
+    UINT32 nl = 0;
+
+    UINT32 nd = 0;
+    UINT32 *arr;
+    int size = 0;
+    UINT32 index;
+    int hammd;
+
+
+    counter->erase();
+    memset(numres, 0, (B+1)*sizeof(*numres));
+
+    split(chunks, Query, m, mplus, b);
+    
+    /* the growing search radius per substring */
+    int s;
+
+    /* current b: for the first mplus substrings it is b, for the rest it is (b-1) */
+    int curb = b;
+
+    for (s = 0; s <= d && n < maxres; s++)
+    {
+        for (int k=0; k<m; k++) {
+            if (k < mplus)
+            curb = b;
+            else
+            curb = b-1;
+            UINT64 chunksk = chunks[k];
+            /* number of bit-strings with s number of 1s */
+            nl += xornum[s+1] - xornum[s];
+
+            /* the bit-string with s number of 1s */
+            UINT64 bitstr = 0;
+            for (int i=0; i<s; i++)
+                /* power[i] stores the location of the i'th 1 */
+                power[i] = i;
+            /* used for stopping criterion (location of (s+1)th 1) */
+            power[s] = curb+1;
+
+            /* bit determines the 1 that should be moving to the left */
+            int bit = s-1;
+
+            /* start from the left-most 1, and move it to the left until
+               it touches another one */
+
+            /* the loop for changing bitstr */
+            while (true)
+            {
+                if (bit != -1)
+                {
+                    bitstr ^= (power[bit] == bit) ? (UINT64)1 << power[bit] : (UINT64)3 << (power[bit]-1);
+                    power[bit]++;
+                    bit--;
+                }
+
+                else
+                { /* bit == -1 */
+                    /* the binary code bitstr is available for processing */
+                    arr = H[k].query(chunksk ^ bitstr, &size); // lookup
+                    if (size)
+                    {   /* the corresponding bucket is not empty */
+                        nd += size;
+                        for (int c = 0; c < size; c++)
+                        {
+                            index = arr[c];
+                            if (!counter->get(index))
+                            { /* if it is not a duplicate */
+                                counter->set(index);
+                                hammd = match(codes.ptr() + (UINT64)index*(B_over_8), Query, B_over_8);
+
+                                nc++;
+                                if (hammd <= D && numres[hammd] < maxres)
+                                   res[hammd * K + numres[hammd]] = index+1;
+
+                                numres[hammd]++;
+                            }
+                    }
+                }
+
+                    /* end of processing */
+                    while (++bit < s && power[bit] == power[bit+1]-1) {
+                    bitstr ^= (UINT64)1 << (power[bit]-1);
+                    power[bit] = bit;
+                    }
+                    if (bit == s)
+                    break;
+                }
+            }
+
+            n = n + numres[s*m+k];
+            if (n >= maxres)
+            break;
+        }
+    }
+        
+    n = 0;
+    for (s = 0; s <= D && n < K; s++ )
+    {
+        for (int c = 0; c < numres[s] && n < K; c++)
+            results[n++] = res[s*K + c];
+    }
+
+}
+
+/* constructor 2 */
+Mihasher::Mihasher(int _B, int _m)
+{
+    B = _B;
+    B_over_8 = B/8;
+    m = _m;
+    b = ceil((double)B/m);
+ 
+    /* assuming that B/2 is large enough radius to include
+       all of the k nearest neighbors */
+    D = ceil(B/2.0);
+    d = ceil((double)D/m);
+   
+    /* mplus is the number of chunks with b bits
+       (m-mplus) is the number of chunks with (b-1) bits */
+    mplus = B - m * (b-1);
+
+    xornum = new UINT32 [d+2];
+    xornum[0] = 0;
+    for (int i=0; i<=d; i++)
+    xornum[i+1] = xornum[i] + choose(b, i);
+    
+    H = new SparseHashtable[m];
+
+    /* H[i].init might fail */
+    for (int i=0; i<mplus; i++)
+        H[i].init(b);
+    for (int i=mplus; i<m; i++)
+        H[i].init(b-1);
+}
+
+/* K setter */
+void Mihasher::setK(int _K)
+{
+    K = _K;
+}
+
+/* desctructor */
+Mihasher::~Mihasher()
+{
+    delete[] xornum;
+    delete[] H;
+}
+
+/* populate tables */
+void Mihasher::populate(cv::Mat & _codes, UINT32 _N, int dim1codes)
+{
+    N = _N;
+    codes = _codes;
+    UINT64 * chunks = new UINT64[m];
+
+    UINT8 * pcodes = codes.ptr();
+    for (UINT64 i=0; i<N; i++, pcodes += dim1codes)
+    {
+        split(chunks, pcodes, m, mplus, b);
+
+        for (int k=0; k<m; k++)
+            H[k].insert(chunks[k], i);
+
+        if (i % (int)ceil(N/1000.0) == 0)
+            fflush(stdout);
+    }
+    
+    delete [] chunks;
+}
+

modules/line_descriptor/src/precomp.hpp

@@ -48,14 +48,18 @@
 #include <opencv2/imgproc.hpp>
 #include <opencv2/features2d.hpp>
 #include <opencv2/highgui.hpp>
+#include "opencv2/core.hpp"
 
-#include "opencv2/line_descriptor.hpp"
-
-#include <math.h>
 #include <iostream>
 #include <map>
-#include <utility>
+#include <stdio.h>
-#include <string>
+#include <string.h>
-#include <typeinfo>
+#include <math.h>
+#include <algorithm>
+#include <bitset>
+
+
+
+#include "opencv2/line_descriptor.hpp"
 
 #endif

modules/line_descriptor/src/sparse_hashtable.cpp

@@ -0,0 +1,40 @@
+#include "precomp.hpp"
+
+const int SparseHashtable::MAX_B = 37;
+
+/* constructor */
+SparseHashtable::SparseHashtable()
+{
+    table = NULL;
+    size = 0;
+    b = 0;
+}
+
+/* initializer */
+int SparseHashtable::init(int _b)
+{
+    b = _b;
+    
+    if (b < 5 || b > MAX_B || b > sizeof(UINT64)*8)
+	return 1;
+    
+    size = UINT64_1 << (b-5);	// size = 2 ^ b
+    table = (BucketGroup*) calloc(size, sizeof(BucketGroup));
+
+    return 0;
+}
+
+/* destructor */
+SparseHashtable::~SparseHashtable () {
+    free(table);
+}
+
+/* insert data */
+void SparseHashtable::insert(UINT64 index, UINT32 data) {
+    table[index >> 5].insert((int)(index % 32), data);
+}
+
+/* query data */
+UINT32* SparseHashtable::query(UINT64 index, int *Size) {
+    return table[index >> 5].query((int)(index % 32), Size);
+}