Fix/vlfeat compliance

src/CMakeLists.txt

@@ -23,6 +23,7 @@ CUDA_ADD_LIBRARY(popsift
         popsift/s_desc_grid.cu popsift/s_desc_grid.h
         popsift/s_desc_igrid.cu popsift/s_desc_igrid.h
         popsift/s_desc_notile.cu popsift/s_desc_notile.h
+        popsift/s_desc_vlfeat.cu popsift/s_desc_vlfeat.h
         popsift/s_desc_norm_rs.h
         popsift/s_desc_norm_l2.h
         popsift/s_desc_normalize.h

src/application/main.cpp

@@ -42,6 +42,7 @@ using namespace std;
 static bool print_dev_info  = false;
 static bool print_time_info = false;
 static bool write_as_uchar  = false;
+static bool write_with_ori  = false;
 static bool dont_write      = false;
 static bool pgmread_loading = false;
 static bool float_mode      = false;
@@ -64,13 +65,13 @@ static void parseargs(int argc, char** argv, popsift::Config& config, string& in
         parameters.add_options()
             ("octaves", value<int>(&config.octaves), "Number of octaves")
             ("levels", value<int>(&config.levels), "Number of levels per octave")
-            ("sigma", value<float>()->notifier([&](float f) { config.setSigma(f); }), "Initial sigma value")
+            ("sigma", value<float>()->notifier([&](float f) { config.setSigma(f); }), "Initial sigma value (default 1.6)")
 
-            ("threshold", value<float>()->notifier([&](float f) { config.setThreshold(f); }), "Contrast threshold")
-            ("edge-threshold", value<float>()->notifier([&](float f) { config.setEdgeLimit(f); }), "On-edge threshold")
+            ("threshold", value<float>()->notifier([&](float f) { config.setThreshold(f); }), "Contrast threshold (default 0.04)")
+            ("edge-threshold", value<float>()->notifier([&](float f) { config.setEdgeLimit(f); }), "On-edge threshold (default 10.0)")
             ("edge-limit", value<float>()->notifier([&](float f) { config.setEdgeLimit(f); }), "On-edge threshold")
             ("downsampling", value<float>()->notifier([&](float f) { config.setDownsampling(f); }), "Downscale width and height of input by 2^N")
-            ("initial-blur", value<float>()->notifier([&](float f) {config.setInitialBlur(f); }), "Assume initial blur, subtract when blurring first time");
+            ("initial-blur", value<float>()->notifier([&](float f) {config.setInitialBlur(f); }), "Assume initial blur, subtract when blurring first time (default 0.5)");
     }
     options_description modes("Modes");
     {
@@ -81,9 +82,9 @@ static void parseargs(int argc, char** argv, popsift::Config& config, string& in
         // "Options are: vlfeat, relative, relative-all, opencv, fixed9, fixed15"
         ("desc-mode", value<std::string>()->notifier([&](const std::string& s) { config.setDescMode(s); }),
         "Choice of descriptor extraction modes:\n"
-        "loop, iloop, grid, igrid, notile\n"
+        "loop, iloop, grid, igrid, notile, vlfeat\n"
 	"Default is loop\n"
-        "loop is OpenCV-like horizontal scanning, computing only valid points, grid extracts only useful points but rounds them, iloop uses linear texture and rotated gradiant fetching. igrid is grid with linear interpolation. notile is like igrid but avoids redundant gradiant fetching.")
+        "loop is OpenCV-like horizontal scanning, computing only valid points, grid extracts only useful points but rounds them, iloop uses linear texture and rotated gradiant fetching. igrid is grid with linear interpolation. notile is like igrid but avoids redundant gradiant fetching, vlfeat replicates vlfeat behaviour.")
         ("popsift-mode", bool_switch()->notifier([&](bool b) { if(b) config.setMode(popsift::Config::PopSift); }),
         "During the initial upscale, shift pixels by 1. In extrema refinement, steps up to 0.6, do not reject points when reaching max iterations, "
         "first contrast threshold is .8 * peak thresh. Shift feature coords octave 0 back to original pos.")
@@ -118,6 +119,7 @@ static void parseargs(int argc, char** argv, popsift::Config& config, string& in
         ("print-time-info", bool_switch(&print_time_info)->default_value(false), "A debug output printing image processing time after load()")
         ("write-as-uchar", bool_switch(&write_as_uchar)->default_value(false), "Output descriptors rounded to int.\n"
          "Scaling to sensible ranges is not automatic, should be combined with --norm-multi=9 or similar")
+        ("write-with-ori", bool_switch(&write_with_ori)->default_value(false), "Output points are written with sigma and orientation.\n")
         ("dont-write", bool_switch(&dont_write)->default_value(false), "Suppress descriptor output")
         ("pgmread-loading", bool_switch(&pgmread_loading)->default_value(false), "Use the old image loader instead of LibDevIL")
         ("float-mode", bool_switch(&float_mode)->default_value(false), "Upload image to GPU as float instead of byte")
@@ -254,7 +256,7 @@ void read_job( SiftJob* job, bool really_write )
         nvtxRangePushA( "Writing features to disk" );
 
         std::ofstream of( "output-features.txt" );
-        feature_list->print( of, write_as_uchar );
+        feature_list->print( of, write_as_uchar, write_with_ori );
     }
     delete feature_list;
 

src/popsift/features.cu

@@ -107,16 +107,16 @@ void FeaturesHost::unpin( )
     cudaHostUnregister( _ori );
 }
 
-void FeaturesHost::print( std::ostream& ostr, bool write_as_uchar ) const
+void FeaturesHost::print( std::ostream& ostr, bool write_as_uchar, bool with_orientation ) const
 {
     for( int i=0; i<size(); i++ ) {
-        _ext[i].print( ostr, write_as_uchar );
+        _ext[i].print( ostr, write_as_uchar, with_orientation );
     }
 }
 
 std::ostream& operator<<( std::ostream& ostr, const FeaturesHost& feature )
 {
-    feature.print( ostr, false );
+    feature.print( ostr, false, false );
     return ostr;
 }
 
@@ -304,13 +304,30 @@ void FeaturesDev::match( FeaturesDev* other )
  * Feature
  *************************************************************/
 
-void Feature::print( std::ostream& ostr, bool write_as_uchar ) const
+void Feature::print( std::ostream& ostr, bool write_as_uchar, bool with_orientation ) const
 {
     float sigval =  1.0f / ( sigma * sigma );
 
     for( int ori=0; ori<num_ori; ori++ ) {
-        ostr << xpos << " " << ypos << " "
-             << sigval << " 0 " << sigval << " ";
+        if( with_orientation )
+        {
+            float dom_ori = orientation[ori];
+            // dom_ori = dom_ori / M_PI2 * 360;
+            // if (dom_ori < 0) dom_ori += 360;
+            if (dom_ori < 0) dom_ori += M_PI2;
+
+            ostr << std::setprecision(6)
+                 << xpos << " " << ypos << " "
+                 << sigma << " "
+                 << dom_ori << " ";
+        }
+        else
+        {
+            ostr << std::setprecision(6)
+                 << xpos << " " << ypos << " "
+                 << sigval << " 0 "
+                 << sigval << " ";
+        }
         if( write_as_uchar ) {
             for( int i=0; i<128; i++ ) {
                 ostr << roundf(desc[ori]->features[i]) << " ";
@@ -328,7 +345,7 @@ void Feature::print( std::ostream& ostr, bool write_as_uchar ) const
 
 std::ostream& operator<<( std::ostream& ostr, const Feature& feature )
 {
-    feature.print( ostr, false );
+    feature.print( ostr, false, false );
     return ostr;
 }
 

src/popsift/features.h

@@ -33,7 +33,7 @@ struct Feature
     float       orientation[ORIENTATION_MAX_COUNT];
     Descriptor* desc[ORIENTATION_MAX_COUNT];
 
-    void print( std::ostream& ostr, bool write_as_uchar ) const;
+    void print( std::ostream& ostr, bool write_as_uchar, bool with_orientation ) const;
 };
 
 std::ostream& operator<<( std::ostream& ostr, const Feature& feature );
@@ -91,7 +91,7 @@ class FeaturesHost : public FeaturesBase
     inline Feature*    getFeatures()    { return _ext; }
     inline Descriptor* getDescriptors() { return _ori; }
 
-    void print( std::ostream& ostr, bool write_as_uchar ) const;
+    void print( std::ostream& ostr, bool write_as_uchar, bool with_orientation ) const;
 
 protected:
     friend class Pyramid;

src/popsift/s_desc_norm_l2.h

@@ -58,15 +58,15 @@ void NormalizeL2::normalize( const float* src_desc, float* dst_desc, const bool
     float norm;
 
     if( threadIdx.x == 0 ) {
-        norm = normf( 128, src_desc );
+        norm = rnormf( 128, src_desc ); // 1/sqrt(sum of squares)
     }
     __syncthreads();
     norm = popsift::shuffle( norm, 0 );
 
-    descr.x = min( descr.x, 0.2f*norm );
-    descr.y = min( descr.y, 0.2f*norm );
-    descr.z = min( descr.z, 0.2f*norm );
-    descr.w = min( descr.w, 0.2f*norm );
+    descr.x = min( descr.x*norm, 0.2f );
+    descr.y = min( descr.y*norm, 0.2f );
+    descr.z = min( descr.z*norm, 0.2f );
+    descr.w = min( descr.w*norm, 0.2f );
 
     norm = descr.x * descr.x
          + descr.y * descr.y
@@ -96,14 +96,14 @@ void NormalizeL2::normalize( const float* src_desc, float* dst_desc, const bool
     norm += popsift::shuffle_down( norm,  2 );
     norm += popsift::shuffle_down( norm,  1 );
     if( threadIdx.x == 0 ) {
-        norm = __fsqrt_rn( norm );
+        norm = __frsqrt_rn( norm );
     }
     norm = popsift::shuffle( norm,  0 );
 
-    descr.x = min( descr.x, 0.2f*norm );
-    descr.y = min( descr.y, 0.2f*norm );
-    descr.z = min( descr.z, 0.2f*norm );
-    descr.w = min( descr.w, 0.2f*norm );
+    descr.x = min( descr.x*norm, 0.2f );
+    descr.y = min( descr.y*norm, 0.2f );
+    descr.z = min( descr.z*norm, 0.2f );
+    descr.w = min( descr.w*norm, 0.2f );
 
     norm = descr.x * descr.x
          + descr.y * descr.y

src/popsift/s_desc_normalize.h

@@ -18,17 +18,17 @@ void normalize_histogram( )
     Descriptor* descs            = dbuf.desc;
     const int   num_orientations = dct.ori_total;
 
-    int offset = blockIdx.x * 32 + threadIdx.y;
+    int offset = blockIdx.x * blockDim.y + threadIdx.y;
 
     // all of these threads are useless
-    if( blockIdx.x * 32 >= num_orientations ) return;
+    if( blockIdx.x * blockDim.y >= num_orientations ) return;
 
-    offset = ( offset < num_orientations ) ? offset
-                                           : num_orientations-1;
-    Descriptor* desc = &descs[offset];
+    offset = min( offset, num_orientations-1 );
+
+    float* features = descs[offset].features;
 
     bool ignoreme = ( offset >= num_orientations );
 
-    T::normalize( desc->features, ignoreme );
+    T::normalize( features, ignoreme );
 }