 void
 __launch_bounds__(BX*BY,1) // max threads, min blocks
 resample_k(const struct work work) {
-    const unsigned idst0 = WORK*(sum(threadIdx)+stride(blockIdx,gridDim)*prod(blockDim));
-    for(unsigned idst=idst0;idst<(idst0+WORK);++idst)
-     {
+    const unsigned ox=threadIdx.x+blockIdx.x*BX;
+    const unsigned oy=(threadIdx.y+blockIdx.y*BY)*WORK;
+    const unsigned idst0=ox+oy*BX;
+
+    for(unsigned idst=idst0;idst<(idst0+WORK*BX);idst+=BX) {
         float lambdas[4];
         unsigned itetrad;
         {

        self->ctx=c;
        memcpy(c->dst_shape,shape,sizeof(c->dst_shape));
        CUTRY(cudaMalloc(&c->dst,sizeof(TPixel)*prod(shape,3)));
    } catch(int) {
        return 0;
    }

        tetrahedron(work.tetrads+i,cubeverts,indexes[i]);

    try {
        const unsigned BX=32*4,BY=1,WORK=32,N=prod(work.dst_shape,3);
        dim3 threads(BX,BY),
             blocks(1,N/BX/BY/WORK);
        CUTRY(cudaGetLastError());
        resample_k<BX,BY,WORK><<<blocks,threads>>>(work);
        CUTRY(cudaGetLastError());
    } catch(int) {
        return 0;
