modules/calib3d/src/opencl/stereobm.cl

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 //
  10 //                           License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2010-2012, Institute Of Software Chinese Academy Of Science, all rights reserved.
  14 // Copyright (C) 2010-2012, Advanced Micro Devices, Inc., all rights reserved.
  15 // Third party copyrights are property of their respective owners.
  16 //
  17 // Redistribution and use in source and binary forms, with or without modification,
  18 // are permitted provided that the following conditions are met:
  19 //
  20 //   * Redistribution's of source code must retain the above copyright notice,
  21 //     this list of conditions and the following disclaimer.
  22 //
  23 //   * Redistribution's in binary form must reproduce the above copyright notice,
  24 //     this list of conditions and the following disclaimer in the documentation
  25 //     and/or other materials provided with the distribution.
  26 //
  27 //   * The name of the copyright holders may not be used to endorse or promote products
  28 //     derived from this software without specific prior written permission.
  29 //
  30 // This software is provided by the copyright holders and contributors as is and
  31 // any express or implied warranties, including, but not limited to, the implied
  32 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  33 // In no event shall the Intel Corporation or contributors be liable for any direct,
  34 // indirect, incidental, special, exemplary, or consequential damages
  35 // (including, but not limited to, procurement of substitute goods or services;
  36 // loss of use, data, or profits; or business interruption) however caused
  37 // and on any theory of liability, whether in contract, strict liability,
  38 // or tort (including negligence or otherwise) arising in any way out of
  39 // the use of this software, even if advised of the possibility of such damage.
  40 //
  41 //M*/
  42
  43 //////////////////////////////////////////////////////////////////////////////////////////////////
  44 ////////////////////////////////////////// stereoBM //////////////////////////////////////////////
  45 //////////////////////////////////////////////////////////////////////////////////////////////////
  46
  47 #define MAX_VAL 32767
  48
  49 #ifndef WSZ
  50 #define WSZ     2
  51 #endif
  52
  53 #define WSZ2    (WSZ / 2)
  54
  55 #ifdef DEFINE_KERNEL_STEREOBM
  56
  57 #define DISPARITY_SHIFT     4
  58 #define FILTERED            ((MIN_DISP - 1) << DISPARITY_SHIFT)
  59
  60 void calcDisp(__local short * cost, __global short * disp, int uniquenessRatio,
  61               __local int * bestDisp, __local int * bestCost, int d, int x, int y, int cols, int rows)
  62 {
  63     int best_disp = *bestDisp, best_cost = *bestCost;
  64     barrier(CLK_LOCAL_MEM_FENCE);
  65
  66     short c = cost[0];
  67     int thresh = best_cost + (best_cost * uniquenessRatio / 100);
  68     bool notUniq = ( (c <= thresh) && (d < (best_disp - 1) || d > (best_disp + 1) ) );
  69
  70     if (notUniq)
  71         *bestCost = FILTERED;
  72     barrier(CLK_LOCAL_MEM_FENCE);
  73
  74     if( *bestCost != FILTERED && x < cols - WSZ2 - MIN_DISP && y < rows - WSZ2 && d == best_disp)
  75     {
  76         int d_aprox = 0;
  77         int yp =0, yn = 0;
  78         if ((0 < best_disp) && (best_disp < NUM_DISP - 1))
  79         {
  80             yp = cost[-2 * BLOCK_SIZE_Y];
  81             yn = cost[2 * BLOCK_SIZE_Y];
  82             d_aprox = yp + yn - 2 * c + abs(yp - yn);
  83         }
  84         disp[0] = (short)(((best_disp + MIN_DISP)*256 + (d_aprox != 0 ? (yp - yn) * 256 / d_aprox : 0) + 15) >> 4);
  85     }
  86 }
  87
  88 short calcCostBorder(__global const uchar * leftptr, __global const uchar * rightptr, int x, int y, int nthread,
  89                      short * costbuf, int *h, int cols, int d, short cost)
  90 {
  91     int head = (*h) % WSZ;
  92     __global const uchar * left, * right;
  93     int idx = mad24(y + WSZ2 * (2 * nthread - 1), cols, x + WSZ2 * (1 - 2 * nthread));
  94     left = leftptr + idx;
  95     right = rightptr + (idx - d);
  96
  97     short costdiff = 0;
  98     if (0 == nthread)
  99     {
 100         #pragma unroll
 101         for (int i = 0; i < WSZ; i++)
 102         {
 103             costdiff += abs( left[0] - right[0] );
 104             left += cols;
 105             right += cols;
 106         }
 107     }
 108     else // (1 == nthread)
 109     {
 110         #pragma unroll
 111         for (int i = 0; i < WSZ; i++)
 112         {
 113             costdiff += abs(left[i] - right[i]);
 114         }
 115     }
 116     cost += costdiff - costbuf[head];
 117     costbuf[head] = costdiff;
 118     *h = head + 1;
 119     return cost;
 120 }
 121
 122 short calcCostInside(__global const uchar * leftptr, __global const uchar * rightptr, int x, int y,
 123                      int cols, int d, short cost_up_left, short cost_up, short cost_left)
 124 {
 125     __global const uchar * left, * right;
 126     int idx = mad24(y - WSZ2 - 1, cols, x - WSZ2 - 1);
 127     left = leftptr + idx;
 128     right = rightptr + (idx - d);
 129     int idx2 = WSZ*cols;
 130
 131     uchar corrner1 = abs(left[0] - right[0]),
 132           corrner2 = abs(left[WSZ] - right[WSZ]),
 133           corrner3 = abs(left[idx2] - right[idx2]),
 134           corrner4 = abs(left[idx2 + WSZ] - right[idx2 + WSZ]);
 135
 136     return cost_up + cost_left - cost_up_left + corrner1 -
 137         corrner2 - corrner3 + corrner4;
 138 }
 139
 140 __kernel void stereoBM(__global const uchar * leftptr,
 141                        __global const uchar * rightptr,
 142                        __global uchar * dispptr, int disp_step, int disp_offset,
 143                        int rows, int cols,                                              // rows, cols of left and right images, not disp
 144                        int textureTreshold, int uniquenessRatio)
 145 {
 146     int lz = get_local_id(0);
 147     int gx = get_global_id(1) * BLOCK_SIZE_X;
 148     int gy = get_global_id(2) * BLOCK_SIZE_Y;
 149
 150     int nthread = lz / NUM_DISP;
 151     int disp_idx = lz % NUM_DISP;
 152
 153     __global short * disp;
 154     __global const uchar * left, * right;
 155
 156     __local short costFunc[2 * BLOCK_SIZE_Y * NUM_DISP];
 157
 158     __local short * cost;
 159     __local int best_disp[2];
 160     __local int best_cost[2];
 161     best_cost[nthread] = MAX_VAL;
 162     best_disp[nthread] = -1;
 163     barrier(CLK_LOCAL_MEM_FENCE);
 164
 165     short costbuf[WSZ];
 166     int head = 0;
 167
 168     int shiftX = WSZ2 + NUM_DISP + MIN_DISP - 1;
 169     int shiftY = WSZ2;
 170
 171     int x = gx + shiftX, y = gy + shiftY, lx = 0, ly = 0;
 172
 173     int costIdx = disp_idx * 2 * BLOCK_SIZE_Y + (BLOCK_SIZE_Y - 1);
 174     cost = costFunc + costIdx;
 175
 176     int tempcost = 0;
 177     if (x < cols - WSZ2 - MIN_DISP && y < rows - WSZ2)
 178     {
 179         if (0 == nthread)
 180         {
 181             #pragma unroll
 182             for (int i = 0; i < WSZ; i++)
 183             {
 184                 int idx = mad24(y - WSZ2, cols, x - WSZ2 + i);
 185                 left = leftptr + idx;
 186                 right = rightptr + (idx - disp_idx);
 187                 short costdiff = 0;
 188                 for(int j = 0; j < WSZ; j++)
 189                 {
 190                     costdiff += abs( left[0] - right[0] );
 191                     left += cols;
 192                     right += cols;
 193                 }
 194                 costbuf[i] = costdiff;
 195             }
 196         }
 197         else // (1 == nthread)
 198         {
 199             #pragma unroll
 200             for (int i = 0; i < WSZ; i++)
 201             {
 202                 int idx = mad24(y - WSZ2 + i, cols, x - WSZ2);
 203                 left = leftptr + idx;
 204                 right = rightptr + (idx - disp_idx);
 205                 short costdiff = 0;
 206                 for (int j = 0; j < WSZ; j++)
 207                 {
 208                     costdiff += abs( left[j] - right[j]);
 209                 }
 210                 tempcost += costdiff;
 211                 costbuf[i] = costdiff;
 212             }
 213         }
 214     }
 215     if (nthread == 1)
 216     {
 217         cost[0] = tempcost;
 218         atomic_min(best_cost + 1, tempcost);
 219     }
 220     barrier(CLK_LOCAL_MEM_FENCE);
 221
 222     if (best_cost[1] == tempcost)
 223          atomic_max(best_disp + 1, disp_idx);
 224     barrier(CLK_LOCAL_MEM_FENCE);
 225
 226     int dispIdx = mad24(gy, disp_step, mad24((int)sizeof(short), gx, disp_offset));
 227     disp = (__global short *)(dispptr + dispIdx);
 228     calcDisp(cost, disp, uniquenessRatio, best_disp + 1, best_cost + 1, disp_idx, x, y, cols, rows);
 229     barrier(CLK_LOCAL_MEM_FENCE);
 230
 231     lx = 1 - nthread;
 232     ly = nthread;
 233
 234     for (int i = 0; i < BLOCK_SIZE_Y * BLOCK_SIZE_X / 2; i++)
 235     {
 236         x = (lx < BLOCK_SIZE_X) ? gx + shiftX + lx : cols;
 237         y = (ly < BLOCK_SIZE_Y) ? gy + shiftY + ly : rows;
 238
 239         best_cost[nthread] = MAX_VAL;
 240         best_disp[nthread] = -1;
 241         barrier(CLK_LOCAL_MEM_FENCE);
 242
 243         costIdx = mad24(2 * BLOCK_SIZE_Y, disp_idx, (BLOCK_SIZE_Y - 1 - ly + lx));
 244         if (0 > costIdx)
 245             costIdx = BLOCK_SIZE_Y - 1;
 246         cost = costFunc + costIdx;
 247         if (x < cols - WSZ2 - MIN_DISP && y < rows - WSZ2)
 248         {
 249             tempcost = (ly * (1 - nthread) + lx * nthread == 0) ?
 250                 calcCostBorder(leftptr, rightptr, x, y, nthread, costbuf, &head, cols, disp_idx, cost[2*nthread-1]) :
 251                 calcCostInside(leftptr, rightptr, x, y, cols, disp_idx, cost[0], cost[1], cost[-1]);
 252         }
 253         cost[0] = tempcost;
 254         atomic_min(best_cost + nthread, tempcost);
 255         barrier(CLK_LOCAL_MEM_FENCE);
 256
 257         if (best_cost[nthread] == tempcost)
 258             atomic_max(best_disp + nthread, disp_idx);
 259         barrier(CLK_LOCAL_MEM_FENCE);
 260
 261         dispIdx = mad24(gy + ly, disp_step, mad24((int)sizeof(short), (gx + lx), disp_offset));
 262         disp = (__global short *)(dispptr + dispIdx);
 263         calcDisp(cost, disp, uniquenessRatio, best_disp + nthread, best_cost + nthread, disp_idx, x, y, cols, rows);
 264
 265         barrier(CLK_LOCAL_MEM_FENCE);
 266
 267         if (lx + nthread - 1 == ly)
 268         {
 269             lx = (lx + nthread + 1) * (1 - nthread);
 270             ly = (ly + 1) * nthread;
 271         }
 272         else
 273         {
 274             lx += nthread;
 275             ly = ly - nthread + 1;
 276         }
 277     }
 278 }
 279 #endif //DEFINE_KERNEL_STEREOBM
 280
 281 //////////////////////////////////////////////////////////////////////////////////////////////////
 282 /////////////////////////////////////// Norm Prefiler ////////////////////////////////////////////
 283 //////////////////////////////////////////////////////////////////////////////////////////////////
 284
 285 __kernel void prefilter_norm(__global unsigned char *input, __global unsigned char *output,
 286                                int rows, int cols, int prefilterCap, int scale_g, int scale_s)
 287 {
 288     // prefilterCap in range 1..63, checked in StereoBMImpl::compute
 289
 290     int x = get_global_id(0);
 291     int y = get_global_id(1);
 292
 293     if(x < cols && y < rows)
 294     {
 295         int cov1 =                                   input[   max(y-1, 0)   * cols + x] * 1 +
 296                   input[y * cols + max(x-1,0)] * 1 + input[      y          * cols + x] * 4 + input[y * cols + min(x+1, cols-1)] * 1 +
 297                                                      input[min(y+1, rows-1) * cols + x] * 1;
 298         int cov2 = 0;
 299         for(int i = -WSZ2; i < WSZ2+1; i++)
 300             for(int j = -WSZ2; j < WSZ2+1; j++)
 301                 cov2 += input[clamp(y+i, 0, rows-1) * cols + clamp(x+j, 0, cols-1)];
 302
 303         int res = (cov1*scale_g - cov2*scale_s)>>10;
 304         res = clamp(res, -prefilterCap, prefilterCap) + prefilterCap;
 305         output[y * cols + x] = res;
 306     }
 307 }
 308
 309
 310 //////////////////////////////////////////////////////////////////////////////////////////////////
 311 ////////////////////////////////////// Sobel Prefiler ////////////////////////////////////////////
 312 //////////////////////////////////////////////////////////////////////////////////////////////////
 313
 314 __kernel void prefilter_xsobel(__global unsigned char *input, __global unsigned char *output,
 315                                int rows, int cols, int prefilterCap)
 316 {
 317     // prefilterCap in range 1..63, checked in StereoBMImpl::compute
 318     int x = get_global_id(0);
 319     int y = get_global_id(1);
 320     if(x < cols && y < rows)
 321     {
 322         if (0 < x && !((y == rows-1) & (rows%2==1) ) )
 323         {
 324             int cov = input[ ((y > 0) ? y-1 : y+1)  * cols + (x-1)] * (-1) + input[ ((y > 0) ? y-1 : y+1)  * cols + ((x<cols-1) ? x+1 : x-1)] * (1) +
 325                       input[              (y)       * cols + (x-1)] * (-2) + input[        (y)             * cols + ((x<cols-1) ? x+1 : x-1)] * (2) +
 326                       input[((y<rows-1)?(y+1):(y-1))* cols + (x-1)] * (-1) + input[((y<rows-1)?(y+1):(y-1))* cols + ((x<cols-1) ? x+1 : x-1)] * (1);
 327
 328             cov = clamp(cov, -prefilterCap, prefilterCap) + prefilterCap;
 329             output[y * cols + x] = cov;
 330         }
 331         else
 332             output[y * cols + x] = prefilterCap;
 333     }
 334 }