 
         dst.ptr(idy)[idx] = d_thinning_lut[log];
 #endif
-    }
-}
+        return true;
+    }
+    return false;
+}
+
+__global__
+void thinning( cv::cuda::PtrStepSzb src, cv::cuda::PtrStepSzb dst )
+{
+    const int block_x = blockIdx.x * V7_WIDTH;
+    const int idx     = block_x + threadIdx.x;
+    const int idy     = blockIdx.y;
+
+    thinning_inner( idx, idy, src, dst, true );
+}
+
+void thinning_and_store( cv::cuda::PtrStepSzb src, cv::cuda::PtrStepSzb dst, uint32_t* edgeCounter, uint32_t edgeMax, int2* edgeCoords )
+{
+    const int block_x = blockIdx.x * V7_WIDTH;
+    const int idx     = block_x + threadIdx.x;
+    const int idy     = blockIdx.y;
+
+    bool keep = thinning_inner( idx, idy, src, dst, false );
+
+    uint32_t mask = __balloc( keep );  // bitfield of warps with results
+    uint32_t ct   = __popc( mask );    // horizontal reduce
+    uint32_t leader = __ffs(mask) - 1; // the highest thread id with indicator==true
+    uint32_t write_index;
+    if( threadIdx.x == leader ) {
+        // leader gets warp's offset from global value and increases it
+        write_index = atomicAdd( edgeCounter, ct );
+    }
+    write_index = __shfl( write_index, leader ); // broadcast warp write index to all
+    write_index += __popc( mask & ((1 << threadIdx.x) - 1) ); // find own write index
+
+    if( keep && write_index < edgeMax ) {
+        edgeCoords[write_index] = make_int4( idx, idy );
+    }
+}
+
+__device__
+void updateXY(const float & dx, const float & dy, int & x, int & y,  float & e, int & stpX, int & stpY)
+{
+    float d = dy / dx;
+    float a = d_abs( d );
+    // stpX = ( dx < 0 ) ? -1 : ( dx == 0 ) ? 0 : 1;
+    // stpY = ( dy < 0 ) ? -1 : ( dy == 0 ) ? 0 : 1;
+    // stpX = ( dx < 0 ) ? -1 : 1;
+    // stpY = ( dy < 0 ) ? -1 : 1;
+    stpX = d_sign( dx );
+    stpY = d_sign( dy );
+    e   += a;
+    x   += stpX;
+    if( e >= 0.5 ) {
+        y += stpY;
+        e -= 1.0f;
+    }
+}
+
+__device__
+bool gradiant_descent_inner( int4&                  out_edge_info,
+                             int2*                  d_edgelist,
+                             uint32_t               edgeCount,
+                             cv::cuda::PtrStepSzb   edges,
+                             // int                    direction,
+                             uint32_t               nmax,
+                             cv::cuda::PtrStepSz16s d_dx,
+                             cv::cuda::PtrStepSz16s d_dy,
+                             int32_t                thrGradient )
+{
+    const int offset    = blockIdx.x * 32 + threadIdx.y;
+    const int direction = threadIdx.x == 0 ? 1 : -1;
+
+    if( offset >= edgeCount ) return false;
+
+    const int idx = d_edgelist[offset].x;
+    const int idy = d_edgelist[offset].y;
+    // const int block_x = blockIdx.x * V7_WIDTH;
+    // const int idx     = block_x + threadIdx.x;
+    // const int idy     = blockIdx.y;
+
+    if( outOfBounds( idx, idy, edges ) ) return false; // should never happen
+
+    if( edges.ptr(idy)[idx] == 0 ) return false; // should never happen
+
+    float  e     = 0.0f;
+    float  dx    = direction * d_dx.ptr(idy)[idx];
+    float  dy    = direction * d_dy.ptr(idy)[idx];
+    const float  adx   = d_abs( dx );
+    const float  ady   = d_abs( dy );
+    size_t n     = 0;
+    int    stpX  = 0;
+    int    stpY  = 0;
+    int    x     = idx;
+    int    y     = idy;
+
+    if( ady > adx ) {
+        updateXY(dy,dx,y,x,e,stpY,stpX);
+    } else {
+        updateXY(dx,dy,x,y,e,stpX,stpY);
+    }
+    n += 1;
+    if ( dx*dx+dy*dy > thrGradient ) {
+        const float dxRef = dx;
+        const float dyRef = dy;
+        const float dx2 = d_dx.ptr(idy)[idx];
+        const float dy2 = d_dy.ptr(idy)[idx];
+        const float compdir = dx2*dxRef+dy2*dyRef;
+        // dir = ( compdir < 0 ) ? -1 : 1;
+        dir = d_sign( compdir );
+        dx = dir * dx2;
+        dy = dir * dy2;
+    }
+    if( ady > adx ) {
+        updateXY(dy,dx,y,x,e,stpY,stpX);
+    } else {
+        updateXY(dx,dy,x,y,e,stpX,stpY);
+    }
+    n += 1;
+
+    if( outOfBounds( x, y, edges ) ) return false;
+
+    ret = edges.ptr(y)[x];
+    if( ret ) {
+        out_edge_info = make_int4( idx, idy, x, y );
+        return true;
+    }
+    
+    while( n <= nmax ) {
+        if( ady > adx ) {
+            updateXY(dy,dx,y,x,e,stpY,stpX);
+        } else {
+            updateXY(dx,dy,x,y,e,stpX,stpY);
+        }
+        n += 1;
+
+        if( outOfBounds( x, y, edges ) ) return false;
+
+        ret = edges.ptr(y)[x];
+        if( ret ) {
+            out_edge_info = make_int4( idx, idy, x, y );
+            return true;
+        }
+
+        if( ady > adx ) {
+            if( outOfBounds( x, y - stpY, edges ) ) return false;
+
+            ret = edges.ptr(y-stpY)[x];
+            if( ret ) {
+                out_edge_info = make_int4( idx, idy, x, y-stpY );
+                return true;
+            }
+        } else {
+            if( outOfBounds( x - stpX, y, edges ) ) return false;
+
+            ret = edges.ptr(y)[x-stpX];
+            if( ret ) {
+                out_edge_info = make_int4( idx, idy, x-stpX, y );
+                return true;
+            }
+        }
+    }
+    return false;
+}
+
+__global__
+void gradiant_descent( int2*                  d_edgelist,
+                       uint32_t               edgeCount,
+                       int4*                  d_new_edgelist,
+                       uint32_t*              d_new_edge_counter,
+                       cv::cuda::PtrStepSzb   edges,
+                       uint32_t               nmax,
+                       cv::cuda::PtrStepSz16s d_dx,
+                       cv::cuda::PtrStepSz16s d_dy,
+                       int32_t                thrGradient )
+{
+    int4 out_edge_info;
+    bool keep = gradiant_descent_inner( out_edge_info,
+                                        d_edgelist,
+                                        edgeCount,
+                                        edges,
+                                        nmax,
+                                        d_dx,
+                                        d_dy,
+                                        thrGradient );
+
+    uint32_t mask = __balloc( keep );  // bitfield of warps with results
+    uint32_t ct   = __popc( mask );    // horizontal reduce
+    uint32_t leader = __ffs(mask) - 1; // the highest thread id with indicator==true
+    uint32_t write_index;
+    if( threadIdx.x == leader ) {
+        // leader gets warp's offset from global value and increases it
+        write_index = atomicAdd( d_new_edge_counter, ct );
+    }
+    write_index = __shfl( write_index, leader ); // broadcast warp write index to all
+    write_index += __popc( mask & ((1 << threadIdx.x) - 1) ); // find own write index
+
+    if( keep && write_index < edgeMax ) {
+        d_new_edgelist[write_index] = out_edge_info;
+    }
+}
+
 
 #if 0
 __global__

    edges.ptr(idy)[idx] = log;
}

__global__
void thinning( cv::cuda::PtrStepSzb src, cv::cuda::PtrStepSzb dst, bool first_run )
{
    const int block_x = blockIdx.x * V7_WIDTH;
    const int idx     = block_x + threadIdx.x;
    const int idy     = blockIdx.y;

    uint8_t log = 0;

    if( src.ptr(idy)[idx] == 0 ) return;

    if( idx >= 1 && idy >=1 && idx <= src.cols-2 && idy <= src.rows-2 ) {
        log |= ( src.ptr(idy-1)[idx-1] != 0 ) ? 0x80 : 0;

    int idx     = block_x + threadIdx.x;
    int idy     = blockIdx.y;

    if( outOfBounds( idx, idy, dst ) ) return;

    int16_t dx = src_dx.ptr(idy)[idx];
    int16_t dy = src_dy.ptr(idy)[idx];
