 
 #define V7_WIDTH    32
 
+__device__
+inline
+bool outOfBounds( int x, int y, const cv::cuda::PtrStepSzb& edges )
+{
+    return ( x < 0 || x >= edges.cols || y < 0 || y >= edges.rows );
+}
+
+__device__
+template <typename T>
+inline T d_abs( T value )
+{
+    return ( ( value < 0 ) ? -value : value );
+}
+
+__device__
+template <typename T>
+inline int d_sign( T value )
+{
+    return ( ( value < 0 ) ? -1 : 1 );
+}
+
 template <class DestType>
 __global__
 void filter_gauss_horiz( cv::cuda::PtrStepSzf          src,

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

    if( idx >= dst.cols ) return;
    if( idy >= dst.rows ) return;

    int16_t dx = src_dx.ptr(idy)[idx];
    int16_t dy = src_dy.ptr(idy)[idx];
    dx = (dx < 0 ) ? -dx : dx;
    dy = (dy < 0 ) ? -dy : dy;
    dst.ptr(idy)[idx] = dx + dy;
}



        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__
