         // coalescent
         float* s = get2DBufferAt(slice, slice_p, depthid, pixid);
         *s = sim;
-    }
-}
-
-__global__ void getBest_kernel(float* slice, int slice_p,
-                               // int* bestDptId, int bestDptId_p,
-                               float* bestDpt, int bestDpt_p,
-                               float* bestSim, int bestSim_p,
-                               int slicesAtTime, int ndepths, int t, int npixs,
-                               int wsh, int width, int height, bool doUsePixelsDepths, int nbest, bool useTcOrRcPixSize,
-                               bool subPixel)
-{
-    int pixid = blockIdx.x * blockDim.x + threadIdx.x;
-    if((pixid < slicesAtTime) && (slicesAtTime * t + pixid < npixs))
-    {
-        float gminSim = -1.1f;
-        for(int bestid = 0; bestid < nbest; bestid++)
-        {
-            float minSim = 1.0f;
-            int minDepthId = 0;
-            for(int depthid = 0; depthid < ndepths; depthid++)
-            {
-                // float sim = slice[pixid*slice_p+depthid];
-                float sim = tex2D(sliceTex, depthid, pixid);
-                if((sim < minSim) && (sim > gminSim))
-                {
-                    minSim = sim;
-                    minDepthId = depthid;
-                }
-            }
-
-            int2 pix = tex2D(pixsTex, pixid, t);
-            patch ptcho;
-            computePatch(ptcho, minDepthId, ndepths, pix, pixid, t, doUsePixelsDepths, useTcOrRcPixSize);
-
-            float3 p = ptcho.p;
-            float floatDepth = size(p - sg_s_rC);
-            float outDepth = floatDepth;
-            float outSim = minSim;
-            // float outDepth = -1.0f;
-            // float outSim   = 1.0f;
-            // subpixel refienement
-            if((subPixel == true) && (minDepthId > 0) && (minDepthId < ndepths - 1))
-            {
-                float3 depths;
-                computePatch(ptcho, minDepthId - 1, ndepths, pix, pixid, t, doUsePixelsDepths, useTcOrRcPixSize);
-
-                p = ptcho.p;
-
-                depths.x = size(p - sg_s_rC);
-                depths.y = floatDepth;
-
-                computePatch(ptcho, minDepthId + 1, ndepths, pix, pixid, t, doUsePixelsDepths, useTcOrRcPixSize);
-
-                p = ptcho.p;
-
-                depths.z = size(p - sg_s_rC);
-
-                float3 sims;
-                sims.x = tex2D(sliceTex, minDepthId - 1, pixid);
-                sims.y = minSim;
-                sims.z = tex2D(sliceTex, minDepthId + 1, pixid);
-
-                float refinedDepth = refineDepthSubPixel(depths, sims);
-
-                if(refinedDepth > 0.0f)
-                {
-                    floatDepth = refinedDepth;
-                    outDepth = refinedDepth;
-                }
-            }
-            float* depth = get2DBufferAt(bestDpt, bestDpt_p, slicesAtTime * t + pixid, bestid);
-            float* sim = get2DBufferAt(bestSim, bestSim_p, slicesAtTime * t + pixid, bestid);
-            *depth = outDepth;
-            *sim = outSim;
-
-            gminSim = minSim;
-        }
-    }
-}
-
-__global__ void getBest_fine_kernel(float* slice, int slice_p,
-                                    // int* bestDptId, int bestDptId_p,
-                                    float* bestDpt, int bestDpt_p,
-                                    float* bestSim, int bestSim_p,
-                                    int slicesAtTime, int ndepths, int t, int npixs,
-                                    int wsh, int width, int height)
-{
-    int pixid = blockIdx.x * blockDim.x + threadIdx.x;
-    if((pixid < slicesAtTime) && (slicesAtTime * t + pixid < npixs))
-    {
-        float minSim = 1.0f;
-        int minDepthId = 0;
-
-        for(int depthid = 0; depthid < ndepths; depthid++)
-        {
-            float sim = tex2D(sliceTex, depthid, pixid);
-            if(sim < minSim)
-            {
-                minSim = sim;
-                minDepthId = depthid;
-            }
-        }
-
-        int2 pix = tex2D(pixsTex, pixid, t);
-        float pixSize;
-        float3 p = computeDepthPoint_fine(pixSize, minDepthId, ndepths, pix, pixid, t);
-        float floatDepth = size(p - sg_s_rC);
-        float outDepth = floatDepth;
-        float outSim = minSim;
-        // float outDepth = -1.0f;
-        // float outSim   = 1.0f;
-        // subpixel refienement
-
-        if((minDepthId > 0) && (minDepthId < ndepths - 1))
-        {
-            // subpixel refine by Stereo Matching with Color-Weighted Correlation, Hierarchical BeliefPropagation, and
-            // Occlusion Handling Qingxiong pami08
-
-            // quadratic polynomial interpolation is used to approximate the cost function between three discrete depth
-            // candidates: d, dA, and dB.
-            // TODO: get formula back from paper as it has been lost by encoding.
-            // d is the discrete depth with the minimal cost, d? ? d ? 1, and d? ? d ? 1. The cost function is approximated as
-            // f?x? ? ax2 ? bx ? c
-
-            float simM1 = tex2D(sliceTex, minDepthId - 1, pixid);
-            float simP1 = tex2D(sliceTex, minDepthId + 1, pixid);
-            float sim1 = minSim;
-            simM1 = (simM1 + 1.0f) / 2.0f;
-            simP1 = (simP1 + 1.0f) / 2.0f;
-            sim1 = (sim1 + 1.0f) / 2.0f;
-
-            if((simM1 > sim1) && (simP1 > sim1))
-            {
-                float dispStep = -((simP1 - simM1) / (2.0f * (simP1 + simM1 - 2.0f * sim1)));
-                p = computeDepthPoint_fine(pixSize, minDepthId - 1, ndepths, pix, pixid, t);
-                float floatDepthM1 = size(p - sg_s_rC);
-                p = computeDepthPoint_fine(pixSize, minDepthId + 1, ndepths, pix, pixid, t);
-                float floatDepthP1 = size(p - sg_s_rC);
-                //-1 : floatDepthM1
-                // 0 : floatDepth
-                //+1 : floatDepthP1
-                // linear function fit
-                // f(x)=a*x+b
-                // floatDepthM1=-a+b
-                // floatDepthP1= a+b
-                // a = b - floatDepthM1
-                // floatDepthP1=2*b-floatDepthM1
-                float b = (floatDepthP1 + floatDepthM1) / 2.0f;
-                float a = b - floatDepthM1;
-                floatDepth = a * dispStep + b;
-                outDepth = floatDepth;
-                outSim = minSim;
-            }
-        }
-
-        assert(false); // TODO: get2DBufferAt
-        bestDpt[slicesAtTime * t + pixid] = outDepth;
-        bestSim[slicesAtTime * t + pixid] = outSim;
     }
 }
 

    }
}

__global__ void downscale_mean_smooth_lab_kernel(
    uchar4* texLab, int texLab_p,
    int width, int height, int scale)
{
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if((x < width) && (y < height))
    {
        float4 t = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
        float sum = 0.0f;
        for(int i = 0; i < scale; i++)
        {
            for(int j = 0; j < scale; j++)
            {
                float4 curPix = 255.0f * tex2D(r4tex, (float)(x * scale + j) + 0.5f, (float)(y * scale + i) + 0.5f);
                t = t + curPix;
                sum += 1.0f;
            }
        }

        t.x = t.x / sum;
        t.y = t.y / sum;
        t.z = t.z / sum;
        t.w = t.w / sum;

        uchar4 tu4 = make_uchar4((unsigned char)t.x, (unsigned char)t.y, (unsigned char)t.z, (unsigned char)t.w);

        *get2DBufferAt(texLab, texLab_p, x, y) = tu4;
    }
}

__global__ void ptsStatForRcDepthMap_kernel(float2* out, int out_p,
                                            float3* pts, int pts_p,
                                            int npts, int width, int height,
                                            int maxNPixSize, int wsh, const float gammaC, const float gammaP)
{
    int ptid = blockIdx.x * blockDim.x + threadIdx.x;

    if(ptid < npts)
    {
        float3 p = pts[ptid];
        float depthP = size(sg_s_rC - p);
        float2 pixf;

        getPixelFor3DPointRC(pixf, p);

        int2 pix = make_int2((int)pixf.x, (int)pixf.y);
        float2 depthVarianceGray = make_float2(0.0f, 0.0f);

        if((pix.x > wsh) && (pix.x < width - wsh) && (pix.y > wsh) && (pix.y < height - wsh))
        {
            float varianceGray = 255.0f * tex2D(r4tex, (float)pix.x + 0.5f, (float)pix.y + 0.5f).w;
            float depth = tex2D(depthsTex, pix.x, pix.y);
            depthVarianceGray = make_float2(depth, varianceGray);
        }

        // TODO!
        out[ptid] = depthVarianceGray;
    }
}

__global__ void getSilhoueteMap_kernel(bool* out, int out_p, int step, int width, int height, const uchar4 maskColorLab)
{
    int x = blockIdx.x * blockDim.x + threadIdx.x;

    float jump = (float)(depthid - ((ndepths - 1) / 2));

    return sg_s_rC + rpv * (depth + pixSize * jump);
}

__global__ void locmin_kernel(float* slice, int slice_p, int ndepths, int slicesAtTime,
