Remove Gaussian blur code alternatives that are never used or didn't

work very well.

Remove the config param ScalingMode (always use default).
Remove fixed scaling code.
Remove code to downscale everything directly from input image.
Remove the narrower gauss filter width called "OpenCV mode".
Remove functions to interpolate from first image plane.
Remove specialized version to create very first level from input image.
Remove Gauss filter tables for direct downscaling using absolute tables.
Removed deprecated scaling mode "OpenCV".

OpenCV was buggy when this code was written. It has improved since then.
Also downscaling by interpolation, which could not be called with any parameter, is removed.
Restructure the calling code for the last 2 pyramid building functions

Move host code for normalized source kernel into kernel's file.
Normalized source mode is only used for the input image. It uses the normalization feature of CUDA textures to scale the input image while creating the first octave.
Simplify the solution with absolute sources.
Return to a solution without shuffle and identical code structure for horizontal and vertical Gaussian filtering.
Host functions to call Gaussian filtering from point textures moved in kernels' code file.
Host functions to call Gaussian filtering from interpolated textures moved in kernels' code file.
Simplified and unified code for absolute source interpolated Gaussian filtering.
Use horiz_from_input_image exclusively for octave 0.
Direct downscaling is not only use for the input image.  Note that initial blur is assumed for every input image, even when it is later interpreted as initially unblurred. That does make a difference, but is apparently recommended.
Extrema refinement modes have more intuitive names and are no longer tied to PopSift vs VLFeat.
(except that the command line parameters of the test code retains the old terms so far)

Author: Carsten Griwodz

src/CMakeLists.txt

@@ -10,10 +10,9 @@ add_library(popsift
         popsift/sift_pyramid.cu popsift/sift_pyramid.h
         popsift/sift_octave.cu popsift/sift_octave.h
         popsift/s_pyramid_build.cu
-        popsift/s_pyramid_build_aa.cu popsift/s_pyramid_build_aa.h
-        popsift/s_pyramid_build_ai.cu popsift/s_pyramid_build_ai.h
-        popsift/s_pyramid_build_ra.cu popsift/s_pyramid_build_ra.h
-        popsift/s_pyramid_fixed.cu
+        popsift/s_pyramid_build_aa.cu
+        popsift/s_pyramid_build_ai.cu
+        popsift/s_pyramid_build_ra.cu
         popsift/sift_extremum.h
         popsift/sift_extremum.cu popsift/s_extrema.cu
         popsift/s_orientation.cu

src/application/main.cpp

@@ -71,28 +71,20 @@ static void parseargs(int argc, char** argv, popsift::Config& config, string& in
         ( "gauss-mode", value<std::string>()->notifier([&](const std::string& s) { config.setGaussMode(s); }),
           popsift::Config::getGaussModeUsage() )
         // "Choice of span (1-sided) for Gauss filters. Default is VLFeat-like computation depending on sigma. "
-        // "Options are: vlfeat, relative, relative-all, opencv, fixed9, fixed15"
+        // "Options are: vlfeat, relative, relative-all, opencv"
         ("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"
-	"Default is loop\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.")
-        ("popsift-mode", bool_switch()->notifier([&](bool b) { if(b) config.setMode(popsift::Config::PopSift); }),
+        ("popsift-mode", bool_switch()->notifier([&](bool b) { if(b) config.setMode(popsift::Config::RefineInOctave); }),
         "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.")
-        ("vlfeat-mode", bool_switch()->notifier([&](bool b) { if(b) config.setMode(popsift::Config::VLFeat); }),
+        ("vlfeat-mode", bool_switch()->notifier([&](bool b) { if(b) config.setMode(popsift::Config::RefineInLevel); }),
         "During the initial upscale, shift pixels by 1. That creates a sharper upscaled image. "
         "In extrema refinement, steps up to 0.6, levels remain unchanged, "
         "do not reject points when reaching max iterations, "
         "first contrast threshold is .8 * peak thresh.")
-        ("opencv-mode", bool_switch()->notifier([&](bool b) { if(b) config.setMode(popsift::Config::OpenCV); }),
-        "During the initial upscale, shift pixels by 0.5. "
-        "In extrema refinement, steps up to 0.5, "
-        "reject points when reaching max iterations, "
-        "first contrast threshold is floor(.5 * peak thresh). "
-        "Computed filter width are lower than VLFeat/PopSift")
-        ("direct-scaling", bool_switch()->notifier([&](bool b) { if(b) config.setScalingMode(popsift::Config::ScaleDirect); }),
-         "Direct each octave from upscaled orig instead of blurred level.")
         ("norm-multi", value<int>()->notifier([&](int i) {config.setNormalizationMultiplier(i); }), "Multiply the descriptor by pow(2,<int>).")
         ( "norm-mode", value<std::string>()->notifier([&](const std::string& s) { config.setNormMode(s); }),
           popsift::Config::getNormModeUsage() )

src/application/match.cpp

@@ -70,28 +70,16 @@ static void parseargs(int argc, char** argv, popsift::Config& config, string& lF
     modes.add_options()
         ( "gauss-mode", value<std::string>()->notifier([&](const std::string& s) { config.setGaussMode(s); }),
           popsift::Config::getGaussModeUsage() )
-        ("desc-mode", value<std::string>()->notifier([&](const std::string& s) { config.setDescMode(s); }),
+        ( "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"
-	"Default is loop\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.")
-        ("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.")
-        ("vlfeat-mode", bool_switch()->notifier([&](bool b) { if(b) config.setMode(popsift::Config::VLFeat); }),
-        "During the initial upscale, shift pixels by 1. That creates a sharper upscaled image. "
-        "In extrema refinement, steps up to 0.6, levels remain unchanged, "
-        "do not reject points when reaching max iterations, "
-        "first contrast threshold is .8 * peak thresh.")
-        ("opencv-mode", bool_switch()->notifier([&](bool b) { if(b) config.setMode(popsift::Config::OpenCV); }),
-        "During the initial upscale, shift pixels by 0.5. "
-        "In extrema refinement, steps up to 0.5, "
-        "reject points when reaching max iterations, "
-        "first contrast threshold is floor(.5 * peak thresh). "
-        "Computed filter width are lower than VLFeat/PopSift")
-        ("direct-scaling", bool_switch()->notifier([&](bool b) { if(b) config.setScalingMode(popsift::Config::ScaleDirect); }),
-         "Direct each octave from upscaled orig instead of blurred level.")
-        ("norm-multi", value<int>()->notifier([&](int i) {config.setNormalizationMultiplier(i); }), "Multiply the descriptor by pow(2,<int>).")
+        ( "popsift-mode", bool_switch()->notifier([&](bool b) { if(b) config.setMode(popsift::Config::RefineInOctave); }),
+        "In extrema refinement, it is possible to move extrema within a level but also between the levels of an octave.")
+        ( "vlfeat-mode", bool_switch()->notifier([&](bool b) { if(b) config.setMode(popsift::Config::RefineInLevel); }),
+        "In extrema refinement, it is possible to move extrema within a level only.")
+        ( "norm-multi", value<int>()->notifier([&](int i) {config.setNormalizationMultiplier(i); }), "Multiply the descriptor by pow(2,<int>).")
         ( "norm-mode", value<std::string>()->notifier([&](const std::string& s) { config.setNormMode(s); }),
           popsift::Config::getNormModeUsage() )
         ( "root-sift", bool_switch()->notifier([&](bool b) { if(b) config.setNormMode(popsift::Config::RootSift); }),

src/popsift/gauss_filter.cu

@@ -67,42 +67,7 @@ void print_gauss_filter_symbol( int columns )
 
     printf( "\n"
             "Gauss tables\n"
-            "      level span sigma : center value -> edge value\n"
-            "      absolute filters octave 0 (compute level 0, all other levels directly from level 0)\n");
-
-    for( int lvl=0; lvl<d_gauss.required_filter_stages; lvl++ ) {
-        int span = d_gauss.abs_o0.span[lvl] + d_gauss.abs_o0.span[lvl] - 1;
-
-        printf("      %d %d %2.6f: ", lvl, span, d_gauss.abs_o0.sigma[lvl] );
-        int m = min( d_gauss.abs_o0.span[lvl], columns );
-        for( int x=0; x<m; x++ ) {
-            printf("%0.8f ", d_gauss.abs_o0.filter[lvl*GAUSS_ALIGN+x] );
-        }
-        if( m < d_gauss.abs_o0.span[lvl] )
-            printf("...\n");
-        else
-            printf("\n");
-    }
-    printf( "\n"
-            "      absolute filters other octaves\n"
-            "      (level 0 via downscaling, all other levels directly from level 0)\n");
-
-    for( int lvl=0; lvl<d_gauss.required_filter_stages; lvl++ ) {
-        int span = d_gauss.abs_oN.span[lvl] + d_gauss.abs_oN.span[lvl] - 1;
-
-        printf("      %d %d %2.6f: ", lvl, span, d_gauss.abs_oN.sigma[lvl] );
-        int m = min( d_gauss.abs_oN.span[lvl], columns );
-        for( int x=0; x<m; x++ ) {
-            printf("%0.8f ", d_gauss.abs_oN.filter[lvl*GAUSS_ALIGN+x] );
-        }
-        if( m < d_gauss.abs_oN.span[lvl] )
-            printf("...\n");
-        else
-            printf("\n");
-    }
-    printf("\n");
-
-    printf("    level 0-filters for direct downscaling\n");
+            "    level 0-filters for direct downscaling\n");
 
     for( int lvl=0; lvl<MAX_OCTAVES; lvl++ ) {
         int span = d_gauss.dd.span[lvl] + d_gauss.dd.span[lvl] - 1;
@@ -187,33 +152,6 @@ void init_filter( const Config& conf,
 
     h_gauss.inc.computeBlurTable( &h_gauss );
 
-    /* abs_o0 :
-     * Gauss table to create octave 0 of the absolute filters directly from
-     * input images.
-     */
-    for( int lvl=0; lvl<h_gauss.required_filter_stages; lvl++ ) {
-        const float sigmaS = sigma0 * pow( 2.0f, (float)(lvl)/(float)levels );
-        h_gauss.abs_o0.sigma[lvl]  = sqrt( fabs( sigmaS * sigmaS - initial_blur * initial_blur ) );
-    }
-
-    h_gauss.abs_o0.computeBlurTable( &h_gauss );
-
-    /* abs_oN :
-     * Gauss tables to create levels 1 and above directly from level 0 of every
-     * octave. Could be used on octave 0, but abs_o0 is better.
-     * Level 0 must be created by other means (downscaling from previous octave,
-     * direct downscaling from input image, ...) before using abs_oN.
-     * 
-     */
-    h_gauss.abs_oN.sigma[0] = 0;
-    for( int lvl=1; lvl<h_gauss.required_filter_stages; lvl++ ) {
-        const float sigmaP = sigma0; // level 0 has already reached sigma0 blur
-        const float sigmaS = sigma0 * pow( 2.0f, (float)(lvl)/(float)levels );
-        h_gauss.abs_oN.sigma[lvl] = sqrt( sigmaS * sigmaS - sigmaP * sigmaP );
-    }
-
-    h_gauss.abs_oN.computeBlurTable( &h_gauss );
-
     /* dd :
      * The direct-downscaling kernels make use of the assumption that downscaling
      * from MAX_LEVEL-3 is identical to applying 2*sigma on the identical image
@@ -224,17 +162,13 @@ void init_filter( const Config& conf,
      * octaves, where it is also good for performance.
      * dd is only for creating level 0 of all octave directly from the input image.
      */
-    for( int oct=0; oct<MAX_OCTAVES; oct++ ) {
-        // sigma * 2^i
-        float oct_sigma = scalbnf( sigma0, oct );
 
-        // subtract initial blur
-        float b = sqrt( fabs( oct_sigma * oct_sigma - initial_blur * initial_blur ) );
+    // subtract initial blur
+    const float b = sqrt( fabs( sigma0 * sigma0 - initial_blur * initial_blur ) );
 
-        // sigma / 2^i
-        h_gauss.dd.sigma[oct] = scalbnf( b, -oct );
-        h_gauss.dd.computeBlurTable( &h_gauss );
-    }
+    // sigma / 2^i
+    h_gauss.dd.sigma[0] = b;
+    h_gauss.dd.computeBlurTable( &h_gauss );
 
     cudaError_t err;
     err = cudaMemcpyToSymbol( d_gauss,
@@ -260,8 +194,6 @@ __host__
 void GaussInfo::clearTables( )
 {
     inc            .clearTables();
-    abs_o0         .clearTables();
-    abs_oN         .clearTables();
     dd             .clearTables();
 }
 
@@ -276,20 +208,10 @@ int GaussInfo::getSpan( float sigma ) const
 {
     switch( _span_mode )
     {
-    case Config::VLFeat_Relative_All :
-        // return GaussInfo::vlFeatRelativeSpan( sigma );
-        return GaussInfo::vlFeatSpan( sigma );
-
     case Config::VLFeat_Compute :
         return GaussInfo::vlFeatSpan( sigma );
     case Config::VLFeat_Relative :
         return GaussInfo::vlFeatRelativeSpan( sigma );
-    case Config::OpenCV_Compute :
-        return GaussInfo::openCVSpan( sigma );
-    case Config::Fixed9 :
-        return 5;
-    case Config::Fixed15 :
-        return 8;
     default :
         stringstream ss;
         ss << "ERROR: The mode for computing Gauss filter scan is invalid";
@@ -317,15 +239,6 @@ int GaussInfo::vlFeatRelativeSpan( float sigma )
     return spn;
 }
 
-__host__
-int GaussInfo::openCVSpan( float sigma )
-{
-    int span = int( roundf( 2.0f * 4.0f * sigma + 1.0f ) ) | 1;
-    span >>= 1;
-    span  += 1;
-    return std::min<int>( span, GAUSS_ALIGN - 1 );
-}
-
 template<int LEVELS>
 __host__
 void GaussTable<LEVELS>::clearTables( )

src/popsift/gauss_filter.h

@@ -61,25 +61,11 @@ struct GaussInfo
      */
     GaussTable<GAUSS_LEVELS> inc;
 
-    /* Compute the 1D Gauss tables for all levels of octave 0.
-     * For octave 0, all of these tables derive from the input
-     * image.
+    /* This is the 1D Gauss table for filtering the input image.
+     * The input image is downscaled and blurred with sigma or by
+     * blurring the input image with 2*sigma and downscaling afterwards.
      */
-    GaussTable<GAUSS_LEVELS> abs_o0;
-
-    /* Compute the 1D Gauss tables for all levels of octaves 1 and up.
-     * Level 0 is empty, since it is created by other means.
-     * All other levels blur from level 0, not considering any
-     * initial blur.
-     */
-    GaussTable<GAUSS_LEVELS> abs_oN;
-
-    /* In theory, level 0 of octave 2 contains the same information
-     * whether it is constructed by downscaling and blurring the
-     * input image with sigma or by blurring the input image with 2*sigma
-     * and downscaling afterwards.
-     */
-    GaussTable<MAX_OCTAVES> dd;
+    GaussTable<1> dd;
 
     __host__
     void clearTables( );
@@ -99,9 +85,6 @@ struct GaussInfo
 
     __host__
     static int vlFeatRelativeSpan( float sigma );
-
-    __host__
-    static int openCVSpan( float sigma );
 };
 
 extern __device__ __constant__ GaussInfo d_gauss;