 #include "nearby_point.h"
 #include "tag_cut.h"
 
+#undef NO_ATOMIC
+
 using namespace std;
 
 namespace popart {

    if( i >= gridNSample ) return;
    if( j >= gridNSample ) return;


    NearbyPoint* nPoint;
    identification::CutSignals* signals;

    {
        const int idx = j * gridNSample + i;
        nPoint  = &point_buffer[idx];
        signals = &sig_buffer[idx * vCutsSize];
    }

    nPoint->point = make_float2( center.x - halfWidth + i*stepSize,
                                center.y - halfWidth + j*stepSize );

    popart::geometry::matrix3x3  mInvT;
    mT.invert( mInvT ); // note: returns false if it fails
    mInvT.condition( nPoint->point );

    nPoint->result   = 0.0f;
    nPoint->resSize  = 0;
    nPoint->readable = true;

 */
__global__
void idComputeResult( const int         tagIndex,
                      NearbyPoint*      point_buffer,
                      const CutStruct*  vCuts,
                      const CutSignals* sig_buffer,
                      const int         vCutsSize )
{
  const size_t gridNSample = tagParam.gridNSample;
  const int grid_i   = blockIdx.y;
  const int grid_j   = blockIdx.z;
  const int grid_idx = grid_j * gridNSample + grid_i;
  const CutSignals* allcut_signals = &sig_buffer[grid_idx * vCutsSize];
  NearbyPoint* const nPoint  = &point_buffer[grid_idx];

#ifdef NO_ATOMIC
  const int numPairs = vCutsSize*(vCutsSize-1)/2;
  for( int myPair = threadIdx.y; myPair < numPairs; myPair += 32 ) {
#else
    int myPair = blockIdx.x * 32 + threadIdx.y;
#endif
    int j      = __float2int_rd( 1.0f + __fsqrt_rd(1.0f+8.0f*myPair) ) / 2;
    int i      = myPair - j*(j-1)/2;

