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_MULTI_DENOISING_INVOKER_HPP__ 43 #define __OPENCV_FAST_NLMEANS_MULTI_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 FastNlMeansMultiDenoisingInvoker : 55 ParallelLoopBody 56 { 57 public: 58 FastNlMeansMultiDenoisingInvoker(const std::vector<Mat>& srcImgs, int imgToDenoiseIndex, 59 int temporalWindowSize, Mat& dst, int template_window_size, 60 int search_window_size, const float *h); 61 62 void operator() (const Range& range) const; 63 64 private: 65 void operator= (const FastNlMeansMultiDenoisingInvoker&); 66 67 int rows_; 68 int cols_; 69 70 Mat& dst_; 71 72 std::vector<Mat> extended_srcs_; 73 Mat main_extended_src_; 74 int border_size_; 75 76 int template_window_size_; 77 int search_window_size_; 78 int temporal_window_size_; 79 80 int template_window_half_size_; 81 int search_window_half_size_; 82 int temporal_window_half_size_; 83 84 typename pixelInfo<WT>::sampleType fixed_point_mult_; 85 int almost_template_window_size_sq_bin_shift; 86 std::vector<WT> almost_dist2weight; 87 88 void calcDistSumsForFirstElementInRow(int i, Array3d<int>& dist_sums, 89 Array4d<int>& col_dist_sums, 90 Array4d<int>& up_col_dist_sums) const; 91 92 void calcDistSumsForElementInFirstRow(int i, int j, int first_col_num, 93 Array3d<int>& dist_sums, Array4d<int>& col_dist_sums, 94 Array4d<int>& up_col_dist_sums) const; 95 }; 96 97 template <typename T, typename IT, typename UIT, typename D, typename WT> 98 FastNlMeansMultiDenoisingInvoker<T, IT, UIT, D, WT>::FastNlMeansMultiDenoisingInvoker( 99 const std::vector<Mat>& srcImgs, 100 int imgToDenoiseIndex, 101 int temporalWindowSize, 102 cv::Mat& dst, 103 int template_window_size, 104 int search_window_size, 105 const float *h) : 106 dst_(dst), extended_srcs_(srcImgs.size()) 107 { 108 CV_Assert(srcImgs.size() > 0); 109 CV_Assert(srcImgs[0].channels() == pixelInfo<T>::channels); 110 111 rows_ = srcImgs[0].rows; 112 cols_ = srcImgs[0].cols; 113 114 template_window_half_size_ = template_window_size / 2; 115 search_window_half_size_ = search_window_size / 2; 116 temporal_window_half_size_ = temporalWindowSize / 2; 117 118 template_window_size_ = template_window_half_size_ * 2 + 1; 119 search_window_size_ = search_window_half_size_ * 2 + 1; 120 temporal_window_size_ = temporal_window_half_size_ * 2 + 1; 121 122 border_size_ = search_window_half_size_ + template_window_half_size_; 123 for (int i = 0; i < temporal_window_size_; i++) 124 copyMakeBorder(srcImgs[imgToDenoiseIndex - temporal_window_half_size_ + i], extended_srcs_[i], 125 border_size_, border_size_, border_size_, border_size_, cv::BORDER_DEFAULT); 126 127 main_extended_src_ = extended_srcs_[temporal_window_half_size_]; 128 const IT max_estimate_sum_value = 129 (IT)temporal_window_size_ * (IT)search_window_size_ * (IT)search_window_size_ * (IT)pixelInfo<T>::sampleMax(); 130 fixed_point_mult_ = (int)std::min<IT>(std::numeric_limits<IT>::max() / max_estimate_sum_value, 131 pixelInfo<WT>::sampleMax()); 132 133 // precalc weight for every possible l2 dist between blocks 134 // additional optimization of precalced weights to replace division(averaging) by binary shift 135 int template_window_size_sq = template_window_size_ * template_window_size_; 136 almost_template_window_size_sq_bin_shift = 0; 137 while (1 << almost_template_window_size_sq_bin_shift < template_window_size_sq) 138 almost_template_window_size_sq_bin_shift++; 139 140 int almost_template_window_size_sq = 1 << almost_template_window_size_sq_bin_shift; 141 double almost_dist2actual_dist_multiplier = (double) almost_template_window_size_sq / template_window_size_sq; 142 143 int max_dist = D::template maxDist<T>(); 144 int almost_max_dist = (int)(max_dist / almost_dist2actual_dist_multiplier + 1); 145 almost_dist2weight.resize(almost_max_dist); 146 147 for (int almost_dist = 0; almost_dist < almost_max_dist; almost_dist++) 148 { 149 double dist = almost_dist * almost_dist2actual_dist_multiplier; 150 almost_dist2weight[almost_dist] = 151 D::template calcWeight<T, WT>(dist, h, fixed_point_mult_); 152 } 153 154 // additional optimization init end 155 if (dst_.empty()) 156 dst_ = Mat::zeros(srcImgs[0].size(), srcImgs[0].type()); 157 } 158 159 template <typename T, typename IT, typename UIT, typename D, typename WT> 160 void FastNlMeansMultiDenoisingInvoker<T, IT, UIT, D, WT>::operator() (const Range& range) const 161 { 162 int row_from = range.start; 163 int row_to = range.end - 1; 164 165 Array3d<int> dist_sums(temporal_window_size_, search_window_size_, search_window_size_); 166 167 // for lazy calc optimization 168 Array4d<int> col_dist_sums(template_window_size_, temporal_window_size_, search_window_size_, search_window_size_); 169 170 int first_col_num = -1; 171 Array4d<int> up_col_dist_sums(cols_, temporal_window_size_, search_window_size_, search_window_size_); 172 173 for (int i = row_from; i <= row_to; i++) 174 { 175 for (int j = 0; j < cols_; j++) 176 { 177 int search_window_y = i - search_window_half_size_; 178 int search_window_x = j - search_window_half_size_; 179 180 // calc dist_sums 181 if (j == 0) 182 { 183 calcDistSumsForFirstElementInRow(i, dist_sums, col_dist_sums, up_col_dist_sums); 184 first_col_num = 0; 185 } 186 else 187 { 188 // calc cur dist_sums using previous dist_sums 189 if (i == row_from) 190 { 191 calcDistSumsForElementInFirstRow(i, j, first_col_num, 192 dist_sums, col_dist_sums, up_col_dist_sums); 193 194 } 195 else 196 { 197 int ay = border_size_ + i; 198 int ax = border_size_ + j + template_window_half_size_; 199 200 int start_by = 201 border_size_ + i - search_window_half_size_; 202 203 int start_bx = 204 border_size_ + j - search_window_half_size_ + template_window_half_size_; 205 206 T a_up = main_extended_src_.at<T>(ay - template_window_half_size_ - 1, ax); 207 T a_down = main_extended_src_.at<T>(ay + template_window_half_size_, ax); 208 209 // copy class member to local variable for optimization 210 int search_window_size = search_window_size_; 211 212 for (int d = 0; d < temporal_window_size_; d++) 213 { 214 Mat cur_extended_src = extended_srcs_[d]; 215 Array2d<int> cur_dist_sums = dist_sums[d]; 216 Array2d<int> cur_col_dist_sums = col_dist_sums[first_col_num][d]; 217 Array2d<int> cur_up_col_dist_sums = up_col_dist_sums[j][d]; 218 for (int y = 0; y < search_window_size; y++) 219 { 220 int* dist_sums_row = cur_dist_sums.row_ptr(y); 221 222 int* col_dist_sums_row = cur_col_dist_sums.row_ptr(y); 223 int* up_col_dist_sums_row = cur_up_col_dist_sums.row_ptr(y); 224 225 const T* b_up_ptr = cur_extended_src.ptr<T>(start_by - template_window_half_size_ - 1 + y); 226 const T* b_down_ptr = cur_extended_src.ptr<T>(start_by + template_window_half_size_ + y); 227 228 for (int x = 0; x < search_window_size; x++) 229 { 230 dist_sums_row[x] -= col_dist_sums_row[x]; 231 232 col_dist_sums_row[x] = up_col_dist_sums_row[x] + 233 D::template calcUpDownDist<T>(a_up, a_down, b_up_ptr[start_bx + x], b_down_ptr[start_bx + x]); 234 235 dist_sums_row[x] += col_dist_sums_row[x]; 236 up_col_dist_sums_row[x] = col_dist_sums_row[x]; 237 } 238 } 239 } 240 } 241 242 first_col_num = (first_col_num + 1) % template_window_size_; 243 } 244 245 // calc weights 246 IT estimation[pixelInfo<T>::channels], weights_sum[pixelInfo<WT>::channels]; 247 for (size_t channel_num = 0; channel_num < pixelInfo<T>::channels; channel_num++) 248 estimation[channel_num] = 0; 249 for (size_t channel_num = 0; channel_num < pixelInfo<WT>::channels; channel_num++) 250 weights_sum[channel_num] = 0; 251 252 for (int d = 0; d < temporal_window_size_; d++) 253 { 254 const Mat& esrc_d = extended_srcs_[d]; 255 for (int y = 0; y < search_window_size_; y++) 256 { 257 const T* cur_row_ptr = esrc_d.ptr<T>(border_size_ + search_window_y + y); 258 259 int* dist_sums_row = dist_sums.row_ptr(d, y); 260 261 for (int x = 0; x < search_window_size_; x++) 262 { 263 int almostAvgDist = dist_sums_row[x] >> almost_template_window_size_sq_bin_shift; 264 265 WT weight = almost_dist2weight[almostAvgDist]; 266 T p = cur_row_ptr[border_size_ + search_window_x + x]; 267 incWithWeight<T, IT, WT>(estimation, weights_sum, weight, p); 268 } 269 } 270 } 271 272 divByWeightsSum<IT, UIT, pixelInfo<T>::channels, pixelInfo<WT>::channels>(estimation, 273 weights_sum); 274 dst_.at<T>(i,j) = saturateCastFromArray<T, IT>(estimation); 275 } 276 } 277 } 278 279 template <typename T, typename IT, typename UIT, typename D, typename WT> 280 inline void FastNlMeansMultiDenoisingInvoker<T, IT, UIT, D, WT>::calcDistSumsForFirstElementInRow( 281 int i, Array3d<int>& dist_sums, Array4d<int>& col_dist_sums, Array4d<int>& up_col_dist_sums) const 282 { 283 int j = 0; 284 285 for (int d = 0; d < temporal_window_size_; d++) 286 { 287 Mat cur_extended_src = extended_srcs_[d]; 288 for (int y = 0; y < search_window_size_; y++) 289 for (int x = 0; x < search_window_size_; x++) 290 { 291 dist_sums[d][y][x] = 0; 292 for (int tx = 0; tx < template_window_size_; tx++) 293 col_dist_sums[tx][d][y][x] = 0; 294 295 int start_y = i + y - search_window_half_size_; 296 int start_x = j + x - search_window_half_size_; 297 298 int* dist_sums_ptr = &dist_sums[d][y][x]; 299 int* col_dist_sums_ptr = &col_dist_sums[0][d][y][x]; 300 int col_dist_sums_step = col_dist_sums.step_size(0); 301 for (int tx = -template_window_half_size_; tx <= template_window_half_size_; tx++) 302 { 303 for (int ty = -template_window_half_size_; ty <= template_window_half_size_; ty++) 304 { 305 int dist = D::template calcDist<T>( 306 main_extended_src_.at<T>(border_size_ + i + ty, border_size_ + j + tx), 307 cur_extended_src.at<T>(border_size_ + start_y + ty, border_size_ + start_x + tx)); 308 309 *dist_sums_ptr += dist; 310 *col_dist_sums_ptr += dist; 311 } 312 col_dist_sums_ptr += col_dist_sums_step; 313 } 314 315 up_col_dist_sums[j][d][y][x] = col_dist_sums[template_window_size_ - 1][d][y][x]; 316 } 317 } 318 } 319 320 template <typename T, typename IT, typename UIT, typename D, typename WT> 321 inline void FastNlMeansMultiDenoisingInvoker<T, IT, UIT, D, WT>::calcDistSumsForElementInFirstRow( 322 int i, int j, int first_col_num, Array3d<int>& dist_sums, 323 Array4d<int>& col_dist_sums, Array4d<int>& up_col_dist_sums) const 324 { 325 int ay = border_size_ + i; 326 int ax = border_size_ + j + template_window_half_size_; 327 328 int start_by = border_size_ + i - search_window_half_size_; 329 int start_bx = border_size_ + j - search_window_half_size_ + template_window_half_size_; 330 331 int new_last_col_num = first_col_num; 332 333 for (int d = 0; d < temporal_window_size_; d++) 334 { 335 Mat cur_extended_src = extended_srcs_[d]; 336 for (int y = 0; y < search_window_size_; y++) 337 for (int x = 0; x < search_window_size_; x++) 338 { 339 dist_sums[d][y][x] -= col_dist_sums[first_col_num][d][y][x]; 340 341 col_dist_sums[new_last_col_num][d][y][x] = 0; 342 int by = start_by + y; 343 int bx = start_bx + x; 344 345 int* col_dist_sums_ptr = &col_dist_sums[new_last_col_num][d][y][x]; 346 for (int ty = -template_window_half_size_; ty <= template_window_half_size_; ty++) 347 { 348 *col_dist_sums_ptr += D::template calcDist<T>( 349 main_extended_src_.at<T>(ay + ty, ax), 350 cur_extended_src.at<T>(by + ty, bx)); 351 } 352 353 dist_sums[d][y][x] += col_dist_sums[new_last_col_num][d][y][x]; 354 355 up_col_dist_sums[j][d][y][x] = col_dist_sums[new_last_col_num][d][y][x]; 356 } 357 } 358 } 359 360 #endif 361
