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 // Intel License Agreement 11 // For Open Source Computer Vision Library 12 // 13 // Copyright (C) 2000, Intel Corporation, all rights reserved. 14 // Third party copyrights are property of their respective icvers. 15 // 16 // Redistribution and use in source and binary forms, with or without modification, 17 // are permitted provided that the following conditions are met: 18 // 19 // * Redistribution's of source code must retain the above copyright notice, 20 // this list of conditions and the following disclaimer. 21 // 22 // * Redistribution's in binary form must reproduce the above copyright notice, 23 // this list of conditions and the following disclaimer in the documentation 24 // and/or other materials provided with the distribution. 25 // 26 // * The name of Intel Corporation may not be used to endorse or promote products 27 // derived from this software without specific prior written permission. 28 // 29 // This software is provided by the copyright holders and contributors "as is" and 30 // any express or implied warranties, including, but not limited to, the implied 31 // warranties of merchantability and fitness for a particular purpose are disclaimed. 32 // In no event shall the Intel Corporation or contributors be liable for any direct, 33 // indirect, incidental, special, exemplary, or consequential damages 34 // (including, but not limited to, procurement of substitute goods or services; 35 // loss of use, data, or profits; or business interruption) however caused 36 // and on any theory of liability, whether in contract, strict liability, 37 // or tort (including negligence or otherwise) arising in any way out of 38 // the use of this software, even if advised of the possibility of such damage. 39 // 40 //M*/ 41 42 #ifndef __OPENCV_FAST_NLMEANS_DENOISING_INVOKER_HPP__ 43 #define __OPENCV_FAST_NLMEANS_DENOISING_INVOKER_HPP__ 44 45 #include "precomp.hpp" 46 #include <limits> 47 48 #include "fast_nlmeans_denoising_invoker_commons.hpp" 49 #include "arrays.hpp" 50 51 using namespace cv; 52 53 template <typename T, typename IT, typename UIT, typename D, typename WT> 54 struct FastNlMeansDenoisingInvoker : 55 public ParallelLoopBody 56 { 57 public: 58 FastNlMeansDenoisingInvoker(const Mat& src, Mat& dst, 59 int template_window_size, int search_window_size, const float *h); 60 61 void operator() (const Range& range) const; 62 63 private: 64 void operator= (const FastNlMeansDenoisingInvoker&); 65 66 const Mat& src_; 67 Mat& dst_; 68 69 Mat extended_src_; 70 int border_size_; 71 72 int template_window_size_; 73 int search_window_size_; 74 75 int template_window_half_size_; 76 int search_window_half_size_; 77 78 typename pixelInfo<WT>::sampleType fixed_point_mult_; 79 int almost_template_window_size_sq_bin_shift_; 80 std::vector<WT> almost_dist2weight_; 81 82 void calcDistSumsForFirstElementInRow( 83 int i, Array2d<int>& dist_sums, 84 Array3d<int>& col_dist_sums, 85 Array3d<int>& up_col_dist_sums) const; 86 87 void calcDistSumsForElementInFirstRow( 88 int i, int j, int first_col_num, 89 Array2d<int>& dist_sums, 90 Array3d<int>& col_dist_sums, 91 Array3d<int>& up_col_dist_sums) const; 92 }; 93 94 inline int getNearestPowerOf2(int value) 95 { 96 int p = 0; 97 while( 1 << p < value) 98 ++p; 99 return p; 100 } 101 102 template <typename T, typename IT, typename UIT, typename D, typename WT> 103 FastNlMeansDenoisingInvoker<T, IT, UIT, D, WT>::FastNlMeansDenoisingInvoker( 104 const Mat& src, Mat& dst, 105 int template_window_size, 106 int search_window_size, 107 const float *h) : 108 src_(src), dst_(dst) 109 { 110 CV_Assert(src.channels() == pixelInfo<T>::channels); 111 112 template_window_half_size_ = template_window_size / 2; 113 search_window_half_size_ = search_window_size / 2; 114 template_window_size_ = template_window_half_size_ * 2 + 1; 115 search_window_size_ = search_window_half_size_ * 2 + 1; 116 117 border_size_ = search_window_half_size_ + template_window_half_size_; 118 copyMakeBorder(src_, extended_src_, border_size_, border_size_, border_size_, border_size_, BORDER_DEFAULT); 119 120 const IT max_estimate_sum_value = 121 (IT)search_window_size_ * (IT)search_window_size_ * (IT)pixelInfo<T>::sampleMax(); 122 fixed_point_mult_ = (int)std::min<IT>(std::numeric_limits<IT>::max() / max_estimate_sum_value, 123 pixelInfo<WT>::sampleMax()); 124 125 // precalc weight for every possible l2 dist between blocks 126 // additional optimization of precalced weights to replace division(averaging) by binary shift 127 CV_Assert(template_window_size_ <= 46340); // sqrt(INT_MAX) 128 int template_window_size_sq = template_window_size_ * template_window_size_; 129 almost_template_window_size_sq_bin_shift_ = getNearestPowerOf2(template_window_size_sq); 130 double almost_dist2actual_dist_multiplier = ((double)(1 << almost_template_window_size_sq_bin_shift_)) / template_window_size_sq; 131 132 int max_dist = D::template maxDist<T>(); 133 int almost_max_dist = (int)(max_dist / almost_dist2actual_dist_multiplier + 1); 134 almost_dist2weight_.resize(almost_max_dist); 135 136 for (int almost_dist = 0; almost_dist < almost_max_dist; almost_dist++) 137 { 138 double dist = almost_dist * almost_dist2actual_dist_multiplier; 139 almost_dist2weight_[almost_dist] = 140 D::template calcWeight<T, WT>(dist, h, fixed_point_mult_); 141 } 142 143 // additional optimization init end 144 if (dst_.empty()) 145 dst_ = Mat::zeros(src_.size(), src_.type()); 146 } 147 148 template <typename T, typename IT, typename UIT, typename D, typename WT> 149 void FastNlMeansDenoisingInvoker<T, IT, UIT, D, WT>::operator() (const Range& range) const 150 { 151 int row_from = range.start; 152 int row_to = range.end - 1; 153 154 // sums of cols anf rows for current pixel p 155 Array2d<int> dist_sums(search_window_size_, search_window_size_); 156 157 // for lazy calc optimization (sum of cols for current pixel) 158 Array3d<int> col_dist_sums(template_window_size_, search_window_size_, search_window_size_); 159 160 int first_col_num = -1; 161 // last elements of column sum (for each element in row) 162 Array3d<int> up_col_dist_sums(src_.cols, search_window_size_, search_window_size_); 163 164 for (int i = row_from; i <= row_to; i++) 165 { 166 for (int j = 0; j < src_.cols; j++) 167 { 168 int search_window_y = i - search_window_half_size_; 169 int search_window_x = j - search_window_half_size_; 170 171 // calc dist_sums 172 if (j == 0) 173 { 174 calcDistSumsForFirstElementInRow(i, dist_sums, col_dist_sums, up_col_dist_sums); 175 first_col_num = 0; 176 } 177 else 178 { 179 // calc cur dist_sums using previous dist_sums 180 if (i == row_from) 181 { 182 calcDistSumsForElementInFirstRow(i, j, first_col_num, 183 dist_sums, col_dist_sums, up_col_dist_sums); 184 } 185 else 186 { 187 int ay = border_size_ + i; 188 int ax = border_size_ + j + template_window_half_size_; 189 190 int start_by = border_size_ + i - search_window_half_size_; 191 int start_bx = border_size_ + j - search_window_half_size_ + template_window_half_size_; 192 193 T a_up = extended_src_.at<T>(ay - template_window_half_size_ - 1, ax); 194 T a_down = extended_src_.at<T>(ay + template_window_half_size_, ax); 195 196 // copy class member to local variable for optimization 197 int search_window_size = search_window_size_; 198 199 for (int y = 0; y < search_window_size; y++) 200 { 201 int * dist_sums_row = dist_sums.row_ptr(y); 202 int * col_dist_sums_row = col_dist_sums.row_ptr(first_col_num, y); 203 int * up_col_dist_sums_row = up_col_dist_sums.row_ptr(j, y); 204 205 const T * b_up_ptr = extended_src_.ptr<T>(start_by - template_window_half_size_ - 1 + y); 206 const T * b_down_ptr = extended_src_.ptr<T>(start_by + template_window_half_size_ + y); 207 208 for (int x = 0; x < search_window_size; x++) 209 { 210 // remove from current pixel sum column sum with index "first_col_num" 211 dist_sums_row[x] -= col_dist_sums_row[x]; 212 213 int bx = start_bx + x; 214 col_dist_sums_row[x] = up_col_dist_sums_row[x] + D::template calcUpDownDist<T>(a_up, a_down, b_up_ptr[bx], b_down_ptr[bx]); 215 216 dist_sums_row[x] += col_dist_sums_row[x]; 217 up_col_dist_sums_row[x] = col_dist_sums_row[x]; 218 } 219 } 220 } 221 222 first_col_num = (first_col_num + 1) % template_window_size_; 223 } 224 225 // calc weights 226 IT estimation[pixelInfo<T>::channels], weights_sum[pixelInfo<WT>::channels]; 227 for (int channel_num = 0; channel_num < pixelInfo<T>::channels; channel_num++) 228 estimation[channel_num] = 0; 229 for (int channel_num = 0; channel_num < pixelInfo<WT>::channels; channel_num++) 230 weights_sum[channel_num] = 0; 231 232 for (int y = 0; y < search_window_size_; y++) 233 { 234 const T* cur_row_ptr = extended_src_.ptr<T>(border_size_ + search_window_y + y); 235 int* dist_sums_row = dist_sums.row_ptr(y); 236 for (int x = 0; x < search_window_size_; x++) 237 { 238 int almostAvgDist = dist_sums_row[x] >> almost_template_window_size_sq_bin_shift_; 239 WT weight = almost_dist2weight_[almostAvgDist]; 240 T p = cur_row_ptr[border_size_ + search_window_x + x]; 241 incWithWeight<T, IT, WT>(estimation, weights_sum, weight, p); 242 } 243 } 244 245 divByWeightsSum<IT, UIT, pixelInfo<T>::channels, pixelInfo<WT>::channels>(estimation, 246 weights_sum); 247 dst_.at<T>(i,j) = saturateCastFromArray<T, IT>(estimation); 248 } 249 } 250 } 251 252 template <typename T, typename IT, typename UIT, typename D, typename WT> 253 inline void FastNlMeansDenoisingInvoker<T, IT, UIT, D, WT>::calcDistSumsForFirstElementInRow( 254 int i, 255 Array2d<int>& dist_sums, 256 Array3d<int>& col_dist_sums, 257 Array3d<int>& up_col_dist_sums) const 258 { 259 int j = 0; 260 261 for (int y = 0; y < search_window_size_; y++) 262 for (int x = 0; x < search_window_size_; x++) 263 { 264 dist_sums[y][x] = 0; 265 for (int tx = 0; tx < template_window_size_; tx++) 266 col_dist_sums[tx][y][x] = 0; 267 268 int start_y = i + y - search_window_half_size_; 269 int start_x = j + x - search_window_half_size_; 270 271 for (int ty = -template_window_half_size_; ty <= template_window_half_size_; ty++) 272 for (int tx = -template_window_half_size_; tx <= template_window_half_size_; tx++) 273 { 274 int dist = D::template calcDist<T>(extended_src_, 275 border_size_ + i + ty, border_size_ + j + tx, 276 border_size_ + start_y + ty, border_size_ + start_x + tx); 277 278 dist_sums[y][x] += dist; 279 col_dist_sums[tx + template_window_half_size_][y][x] += dist; 280 } 281 282 up_col_dist_sums[j][y][x] = col_dist_sums[template_window_size_ - 1][y][x]; 283 } 284 } 285 286 template <typename T, typename IT, typename UIT, typename D, typename WT> 287 inline void FastNlMeansDenoisingInvoker<T, IT, UIT, D, WT>::calcDistSumsForElementInFirstRow( 288 int i, int j, int first_col_num, 289 Array2d<int>& dist_sums, 290 Array3d<int>& col_dist_sums, 291 Array3d<int>& up_col_dist_sums) const 292 { 293 int ay = border_size_ + i; 294 int ax = border_size_ + j + template_window_half_size_; 295 296 int start_by = border_size_ + i - search_window_half_size_; 297 int start_bx = border_size_ + j - search_window_half_size_ + template_window_half_size_; 298 299 int new_last_col_num = first_col_num; 300 301 for (int y = 0; y < search_window_size_; y++) 302 for (int x = 0; x < search_window_size_; x++) 303 { 304 dist_sums[y][x] -= col_dist_sums[first_col_num][y][x]; 305 306 col_dist_sums[new_last_col_num][y][x] = 0; 307 int by = start_by + y; 308 int bx = start_bx + x; 309 for (int ty = -template_window_half_size_; ty <= template_window_half_size_; ty++) 310 col_dist_sums[new_last_col_num][y][x] += D::template calcDist<T>(extended_src_, ay + ty, ax, by + ty, bx); 311 312 dist_sums[y][x] += col_dist_sums[new_last_col_num][y][x]; 313 up_col_dist_sums[j][y][x] = col_dist_sums[new_last_col_num][y][x]; 314 } 315 } 316 317 #endif 318
