         /* If the translation of the keypoint is big, move the keypoint
          * and re-iterate the computation. Otherwise we are all set.
          */
-        if( fabs(d.z) < 0.5f && fabs(d.y) < 0.5f && fabs(d.x) < 0.5f) break;
-
-        float3 t;
-
-        t.x = ((d.x >= 0.5f && n.x < width-2) ?  1 : 0 )
-             + ((d.x <= -0.5f && n.x > 1)? -1 : 0 );
-
-        t.y = ((d.y >= 0.5f && n.y < height-2)  ?  1 : 0 )
-             + ((d.y <= -0.5f && n.y > 1) ? -1 : 0 );
-
-        t.z = ((d.z >= 0.5f && n.z < maxlevel-1)  ?  1 : 0 )
-             + ((d.z <= -0.5f && n.z > 1) ? -1 : 0 );
-
-        n.x += t.x;
-        n.y += t.y;
-        n.z += t.z;
+        const bool retval = f.refine( d, n, width, height, maxlevel );
+
+        if( not retval ) break;
     } /* go to next iter */
 
+#ifdef PRINT_EXTREMA_DEBUG_INFO
+    if (iter >= MAX_ITERATIONS) {
+        if( sift_mode == Config::OpenCV ) {
+            printf("Found an extremum at %d %d (o=%d,l=%d) - rejected in refinement, was moved to l:%d (%d,%d)\n", x+1, y+1, debug_octave, level, n.z, n.x, n.y );
+            return false;
+        } else {
+            printf("Found an extremum at %d %d (o=%d,l=%d) - refined to (%d,%d) (o=%d,l=%d)\n", x+1, y+1, debug_octave, level, n.x, n.y, debug_octave, n.z );
+        }
+    }
+#else // PRINT_EXTREMA_DEBUG_INFO
     /* ensure convergence of interpolation */
-    if (iter >= MAX_ITERATIONS) {
-#ifdef PRINT_EXTREMA_DEBUG_INFO
-        printf("Found an extremum at %d %d (o=%d) - rejected in refinement, was moved to l:%d (%d,%d)\n", x+1, y+1, debug_octave, n.z, n.x, n.y );
+    if( sift_mode == Config::OpenCV && iter >= MAX_ITERATIONS) {
+        return false;
+    }
 #endif // PRINT_EXTREMA_DEBUG_INFO
-        return false;
-    }
 
     // float contr   = v + 0.5f * (D.x * d.x + D.y * d.y + D.z * d.z);
     float contr   = v + scalbnf( D.x * d.x + D.y * d.y + D.z * d.z , -1 );

    if( fabs(contr) < scalbnf( d_threshold, 1 ) )
    {
#ifdef PRINT_EXTREMA_DEBUG_INFO
        printf("Found an extremum at %d %d (o=%d) - 2nd peak tresh failed\n", x+1, y+1, debug_octave );
#endif // PRINT_EXTREMA_DEBUG_INFO
        return false;
    }

    return true;
}

template<int sift_mode>
class ModeFunctions
{
public:
    inline __device__
    bool first_contrast_ok( const float val ) const;

    inline __device__
    bool refine( const float3& d, int3& n, const int width, const int height, const int maxlevel );
};

template<>
class ModeFunctions<Config::VLFeat>
{
public:
    inline __device__
    bool first_contrast_ok( const float val ) const
    {
        return ( fabs( val ) >= 0.8 * 2.0 * d_threshold );
    }

    inline __device__
    bool refine( const float3& d, int3& n, const int width, const int height, const int maxlevel )
    {
        float3 t;

        t.x = ((d.x >= 0.6f && n.x < width-2) ?  1 : 0 )
            + ((d.x <= -0.6f && n.x > 1)? -1 : 0 );

        t.y = ((d.y >= 0.6f && n.y < height-2)  ?  1 : 0 )
            + ((d.y <= -0.6f && n.y > 1) ? -1 : 0 );

        // t.z = ((d.z >= 0.6f && n.z < maxlevel-1)  ?  1 : 0 )
            // + ((d.z <= -0.6f && n.z > 1) ? -1 : 0 );

        if( t.x == 0 && t.y == 0 ) return false;

        n.x += t.x;
        n.y += t.y;
        // n.z += t.z; - VLFeat is not changing levels !!!

        return true;
    }
};

template<>
class ModeFunctions<Config::OpenCV>
{
public:
    inline __device__
    bool first_contrast_ok( const float val ) const
    {
        return ( fabs( val ) >= d_threshold );
    }

    inline __device__
    bool refine( const float3& d, int3& n, const int width, const int height, const int maxlevel )
    {
        float3 t;

        t.x = ((d.x >= 0.5f && n.x < width-2) ?  1 : 0 )
            + ((d.x <= -0.5f && n.x > 1)? -1 : 0 );

        t.y = ((d.y >= 0.5f && n.y < height-2)  ?  1 : 0 )
            + ((d.y <= -0.5f && n.y > 1) ? -1 : 0 );

        t.z = ((d.z >= 0.5f && n.z < maxlevel-1)  ?  1 : 0 )
            + ((d.z <= -0.5f && n.z > 1) ? -1 : 0 );

        if( t.x == 0 && t.y == 0 && t.z == 0 ) return false;

        n.x += t.x;
        n.y += t.y;
        n.z += t.z;

        return true;
    }
};

template<>
class ModeFunctions<Config::PopSift>
{
public:
    inline __device__
    bool first_contrast_ok( const float val ) const
    {
        return ( fabs( val ) >= 1.6 * d_threshold );
    }

    inline __device__
    bool refine( float3& d, int3& n, int width, int height, int maxlevel )
    {
        float3 t;

        t.x = ((d.x >= 0.6f && n.x < width-2) ?  1 : 0 )
            + ((d.x <= -0.6f && n.x > 1)? -1 : 0 );

        t.y = ((d.y >= 0.6f && n.y < height-2)  ?  1 : 0 )
            + ((d.y <= -0.6f && n.y > 1) ? -1 : 0 );

        t.z = ((d.z >= 0.6f && n.z < maxlevel-1)  ?  1 : 0 )
            + ((d.z <= -0.6f && n.z > 1) ? -1 : 0 );

        if( t.x == 0 && t.y == 0 && t.z == 0 ) return false;

        n.x += t.x;
        n.y += t.y;
        n.z += t.z;

        return true;
    }
};

template<int sift_mode>
__device__
bool find_extrema_in_dog_v6_sub( cudaTextureObject_t dog,
                                 int                 debug_octave,
