modules/ocl/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 // @Authors
  18 //    Jia Haipeng, jiahaipeng95@gmail.com
  19 //    Sen Liu, swjtuls1987@126.com
  20 //    Peng Xiao, pengxiao@outlook.com
  21 //
  22 // Redistribution and use in source and binary forms, with or without modification,
  23 // are permitted provided that the following conditions are met:
  24 //
  25 //   * Redistribution's of source code must retain the above copyright notice,
  26 //     this list of conditions and the following disclaimer.
  27 //
  28 //   * Redistribution's in binary form must reproduce the above copyright notice,
  29 //     this list of conditions and the following disclaimer in the documentation
  30 //     and/or other oclMaterials provided with the distribution.
  31 //
  32 //   * The name of the copyright holders may not be used to endorse or promote products
  33 //     derived from this software without specific prior written permission.
  34 //
  35 // This software is provided by the copyright holders and contributors as is and
  36 // any express or implied warranties, including, but not limited to, the implied
  37 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  38 // In no event shall the Intel Corporation or contributors be liable for any direct,
  39 // indirect, incidental, special, exemplary, or consequential damages
  40 // (including, but not limited to, procurement of substitute goods or services;
  41 // loss of use, data, or profits; or business interruption) however caused
  42 // and on any theory of liability, whether in contract, strict liability,
  43 // or tort (including negligence or otherwise) arising in any way out of
  44 // the use of this software, even if advised of the possibility of such damage.
  45 //
  46 //M*/
  47
  48 #define ROWSperTHREAD 21     // the number of rows a thread will process
  49 #define BLOCK_W       128    // the thread block width (464)
  50 #define N_DISPARITIES 8
  51
  52 #define STEREO_MIND 0                    // The minimum d range to check
  53 #define STEREO_DISP_STEP N_DISPARITIES   // the d step, must be <= 1 to avoid aliasing
  54
  55 #ifndef radius
  56 #define radius 64
  57 #endif
  58
  59 unsigned int CalcSSD(__local unsigned int *col_ssd)
  60 {
  61     unsigned int cache = col_ssd[0];
  62
  63 #pragma unroll
  64     for(int i = 1; i <= (radius << 1); i++)
  65         cache += col_ssd[i];
  66
  67     return cache;
  68 }
  69
  70 uint2 MinSSD(__local unsigned int *col_ssd)
  71 {
  72     unsigned int ssd[N_DISPARITIES];
  73     const int win_size = (radius << 1);
  74
  75     //See above:  #define COL_SSD_SIZE (BLOCK_W + WIN_SIZE)
  76     ssd[0] = CalcSSD(col_ssd + 0 * (BLOCK_W + win_size));
  77     ssd[1] = CalcSSD(col_ssd + 1 * (BLOCK_W + win_size));
  78     ssd[2] = CalcSSD(col_ssd + 2 * (BLOCK_W + win_size));
  79     ssd[3] = CalcSSD(col_ssd + 3 * (BLOCK_W + win_size));
  80     ssd[4] = CalcSSD(col_ssd + 4 * (BLOCK_W + win_size));
  81     ssd[5] = CalcSSD(col_ssd + 5 * (BLOCK_W + win_size));
  82     ssd[6] = CalcSSD(col_ssd + 6 * (BLOCK_W + win_size));
  83     ssd[7] = CalcSSD(col_ssd + 7 * (BLOCK_W + win_size));
  84
  85     unsigned int mssd = min(min(min(ssd[0], ssd[1]), min(ssd[4], ssd[5])), min(min(ssd[2], ssd[3]), min(ssd[6], ssd[7])));
  86
  87     int bestIdx = 0;
  88
  89     for (int i = 0; i < N_DISPARITIES; i++)
  90     {
  91         if (mssd == ssd[i])
  92             bestIdx = i;
  93     }
  94
  95     return (uint2)(mssd, bestIdx);
  96 }
  97
  98 void StepDown(int idx1, int idx2, __global unsigned char* imageL,
  99               __global unsigned char* imageR, int d,   __local unsigned int *col_ssd)
 100 {
 101     uint8 imgR1 = convert_uint8(vload8(0, imageR + (idx1 - d - 7)));
 102     uint8 imgR2 = convert_uint8(vload8(0, imageR + (idx2 - d - 7)));
 103     uint8 diff1 = (uint8)(imageL[idx1]) - imgR1;
 104     uint8 diff2 = (uint8)(imageL[idx2]) - imgR2;
 105     uint8 res = diff2 * diff2 - diff1 * diff1;
 106     const int win_size = (radius << 1);
 107     col_ssd[0 * (BLOCK_W + win_size)] += res.s7;
 108     col_ssd[1 * (BLOCK_W + win_size)] += res.s6;
 109     col_ssd[2 * (BLOCK_W + win_size)] += res.s5;
 110     col_ssd[3 * (BLOCK_W + win_size)] += res.s4;
 111     col_ssd[4 * (BLOCK_W + win_size)] += res.s3;
 112     col_ssd[5 * (BLOCK_W + win_size)] += res.s2;
 113     col_ssd[6 * (BLOCK_W + win_size)] += res.s1;
 114     col_ssd[7 * (BLOCK_W + win_size)] += res.s0;
 115 }
 116
 117 void InitColSSD(int x_tex, int y_tex, int im_pitch, __global unsigned char* imageL,
 118                 __global unsigned char* imageR, int d,
 119                  __local unsigned int *col_ssd)
 120 {
 121     uint8 leftPixel1;
 122     uint8 diffa = 0;
 123     int idx = y_tex * im_pitch + x_tex;
 124     const int win_size = (radius << 1);
 125     for(int i = 0; i < (win_size + 1); i++)
 126     {
 127         leftPixel1 = (uint8)(imageL[idx]);
 128         uint8 imgR = convert_uint8(vload8(0, imageR + (idx - d - 7)));
 129         uint8 res = leftPixel1 - imgR;
 130         diffa += res * res;
 131
 132         idx += im_pitch;
 133     }
 134     //See above:  #define COL_SSD_SIZE (BLOCK_W + WIN_SIZE)
 135     col_ssd[0 * (BLOCK_W + win_size)] = diffa.s7;
 136     col_ssd[1 * (BLOCK_W + win_size)] = diffa.s6;
 137     col_ssd[2 * (BLOCK_W + win_size)] = diffa.s5;
 138     col_ssd[3 * (BLOCK_W + win_size)] = diffa.s4;
 139     col_ssd[4 * (BLOCK_W + win_size)] = diffa.s3;
 140     col_ssd[5 * (BLOCK_W + win_size)] = diffa.s2;
 141     col_ssd[6 * (BLOCK_W + win_size)] = diffa.s1;
 142     col_ssd[7 * (BLOCK_W + win_size)] = diffa.s0;
 143 }
 144
 145 __kernel void stereoKernel(__global unsigned char *left, __global unsigned char *right,
 146                            __global unsigned int *cminSSDImage, int cminSSD_step,
 147                            __global unsigned char *disp, int disp_step,int cwidth, int cheight,
 148                            int img_step, int maxdisp,
 149                            __local unsigned int *col_ssd_cache)
 150 {
 151     __local unsigned int *col_ssd = col_ssd_cache + get_local_id(0);
 152     __local unsigned int *col_ssd_extra = get_local_id(0) < (radius << 1) ? col_ssd + BLOCK_W : 0;
 153
 154     int X = get_group_id(0) * BLOCK_W + get_local_id(0) + maxdisp + radius;
 155
 156 #define Y (get_group_id(1) * ROWSperTHREAD + radius)
 157
 158     __global unsigned int* minSSDImage = cminSSDImage + X + Y * cminSSD_step;
 159     __global unsigned char* disparImage = disp + X + Y * disp_step;
 160
 161     int end_row = ROWSperTHREAD < (cheight - Y) ? ROWSperTHREAD:(cheight - Y);
 162     int y_tex;
 163     int x_tex = X - radius;
 164
 165     if (x_tex >= cwidth)
 166         return;
 167
 168     for(int d = STEREO_MIND; d < maxdisp; d += STEREO_DISP_STEP)
 169     {
 170         y_tex = Y - radius;
 171
 172         InitColSSD(x_tex, y_tex, img_step, left, right, d, col_ssd);
 173         if (col_ssd_extra > 0)
 174             if (x_tex + BLOCK_W < cwidth)
 175                 InitColSSD(x_tex + BLOCK_W, y_tex, img_step, left, right, d, col_ssd_extra);
 176
 177         barrier(CLK_LOCAL_MEM_FENCE); //before MinSSD function
 178
 179         uint2 minSSD = MinSSD(col_ssd);
 180         if (X < cwidth - radius && Y < cheight - radius)
 181         {
 182             if (minSSD.x < minSSDImage[0])
 183             {
 184                 disparImage[0] = (unsigned char)(d + minSSD.y);
 185                 minSSDImage[0] = minSSD.x;
 186             }
 187         }
 188
 189         for(int row = 1; row < end_row; row++)
 190         {
 191             int idx1 = y_tex * img_step + x_tex;
 192             int idx2 = min(y_tex + ((radius << 1) + 1), cheight - 1) * img_step + x_tex;
 193
 194             barrier(CLK_LOCAL_MEM_FENCE);
 195
 196             StepDown(idx1, idx2, left, right, d, col_ssd);
 197             if (col_ssd_extra > 0)
 198                 if (x_tex + BLOCK_W < cwidth)
 199                     StepDown(idx1, idx2, left + BLOCK_W, right + BLOCK_W, d, col_ssd_extra);
 200
 201             barrier(CLK_LOCAL_MEM_FENCE);
 202
 203             uint2 minSSD = MinSSD(col_ssd);
 204             if (X < cwidth - radius && row < cheight - radius - Y)
 205             {
 206                 int idx = row * cminSSD_step;
 207                 if (minSSD.x < minSSDImage[idx])
 208                 {
 209                     disparImage[disp_step * row] = (unsigned char)(d + minSSD.y);
 210                     minSSDImage[idx] = minSSD.x;
 211                 }
 212             }
 213
 214             y_tex++;
 215         } // for row loop
 216     } // for d loop
 217 }
 218 //////////////////////////////////////////////////////////////////////////////////////////////////
 219 //////////////////////////// Sobel Prefiler (signal channel)//////////////////////////////////////
 220 //////////////////////////////////////////////////////////////////////////////////////////////////
 221
 222 __kernel void prefilter_xsobel(__global unsigned char *input, __global unsigned char *output,
 223                                int rows, int cols, int prefilterCap)
 224 {
 225     int x = get_global_id(0);
 226     int y = get_global_id(1);
 227
 228     if(x < cols && y < rows)
 229     {
 230         int cov = input[(y-1) * cols + (x-1)] * (-1) + input[(y-1) * cols + (x+1)] * (1) +
 231                   input[(y)   * cols + (x-1)] * (-2) + input[(y)   * cols + (x+1)] * (2) +
 232                   input[(y+1) * cols + (x-1)] * (-1) + input[(y+1) * cols + (x+1)] * (1);
 233
 234         cov = min(min(max(-prefilterCap, cov), prefilterCap) + prefilterCap, 255);
 235         output[y * cols + x] = cov & 0xFF;
 236     }
 237 }
 238
 239
 240 //////////////////////////////////////////////////////////////////////////////////////////////////
 241 /////////////////////////////////// Textureness filtering ////////////////////////////////////////
 242 //////////////////////////////////////////////////////////////////////////////////////////////////
 243
 244 float sobel(__global unsigned char *input, int x, int y, int rows, int cols)
 245 {
 246     float conv = 0;
 247     int y1 = y==0? 0 : y-1;
 248     int x1 = x==0? 0 : x-1;
 249     if(x < cols && y < rows && x > 0 && y > 0)
 250     {
 251         conv = (float)input[(y1)  * cols + (x1)] * (-1) + (float)input[(y1)  * cols + (x+1)] * (1) +
 252                (float)input[(y)   * cols + (x1)] * (-2) + (float)input[(y)   * cols + (x+1)] * (2) +
 253                (float)input[(y+1) * cols + (x1)] * (-1) + (float)input[(y+1) * cols + (x+1)] * (1);
 254
 255     }
 256     return fabs(conv);
 257 }
 258
 259 float CalcSums(__local float *cols, __local float *cols_cache, int winsz)
 260 {
 261     float cache = 0;
 262     float cache2 = 0;
 263     int winsz2 = winsz/2;
 264
 265     int x = get_local_id(0);
 266     int group_size_x = get_local_size(0);
 267
 268     for(int i = 1; i <= winsz2; i++)
 269         cache += cols[i];
 270
 271     cols_cache[0] = cache;
 272
 273     barrier(CLK_LOCAL_MEM_FENCE);
 274
 275     if (x < group_size_x - winsz2)
 276         cache2 = cols_cache[winsz2];
 277     else
 278         for(int i = winsz2 + 1; i < winsz; i++)
 279             cache2 += cols[i];
 280
 281     return cols[0] + cache + cache2;
 282 }
 283
 284 #define RpT (2 * ROWSperTHREAD)  // got experimentally
 285 __kernel void textureness_kernel(__global unsigned char *disp, int disp_rows, int disp_cols,
 286                                  int disp_step, __global unsigned char *input, int input_rows,
 287                                  int input_cols,int winsz, float threshold,
 288                                  __local float *cols_cache)
 289 {
 290     int winsz2 = winsz/2;
 291     int n_dirty_pixels = (winsz2) * 2;
 292
 293     int local_id_x = get_local_id(0);
 294     int group_size_x = get_local_size(0);
 295     int group_id_y = get_group_id(1);
 296
 297     __local float *cols = cols_cache + group_size_x + local_id_x;
 298     __local float *cols_extra = local_id_x < n_dirty_pixels ? cols + group_size_x : 0;
 299
 300     int x = get_global_id(0);
 301     int beg_row = group_id_y * RpT;
 302     int end_row = min(beg_row + RpT, disp_rows);
 303
 304 //   if (x < disp_cols)
 305 //   {
 306     int y = beg_row;
 307
 308     float sum = 0;
 309     float sum_extra = 0;
 310
 311     for(int i = y - winsz2; i <= y + winsz2; ++i)
 312     {
 313         sum += sobel(input, x - winsz2, i, input_rows, input_cols);
 314         if (cols_extra)
 315             sum_extra += sobel(input, x + group_size_x - winsz2, i, input_rows, input_cols);
 316     }
 317     *cols = sum;
 318     if (cols_extra)
 319         *cols_extra = sum_extra;
 320
 321     barrier(CLK_LOCAL_MEM_FENCE);
 322
 323     float sum_win = CalcSums(cols, cols_cache + local_id_x, winsz) * 255;
 324     if (sum_win < threshold)
 325         disp[y * disp_step + x] = 0;
 326
 327     barrier(CLK_LOCAL_MEM_FENCE);
 328
 329     for(int y = beg_row + 1; y < end_row; ++y)
 330     {
 331         sum = sum - sobel(input, x - winsz2, y - winsz2 - 1, input_rows, input_cols) +
 332               sobel(input, x - winsz2, y + winsz2, input_rows, input_cols);
 333         *cols = sum;
 334
 335         if (cols_extra)
 336         {
 337             sum_extra = sum_extra - sobel(input, x + group_size_x - winsz2, y - winsz2 - 1,input_rows, input_cols)
 338                         + sobel(input, x + group_size_x - winsz2, y + winsz2, input_rows, input_cols);
 339             *cols_extra = sum_extra;
 340         }
 341
 342         barrier(CLK_LOCAL_MEM_FENCE);
 343         float sum_win = CalcSums(cols, cols_cache + local_id_x, winsz) * 255;
 344         if (sum_win < threshold)
 345             disp[y * disp_step + x] = 0;
 346
 347         barrier(CLK_LOCAL_MEM_FENCE);
 348     }
 349     //  }
 350 }