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>
FastNlMeansMultiDenoisingInvoker(const std::vector<Mat> & srcImgs,int imgToDenoiseIndex,int temporalWindowSize,cv::Mat & dst,int template_window_size,int search_window_size,const float * h)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>
operator ()(const Range & range) const160 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>
calcDistSumsForFirstElementInRow(int i,Array3d<int> & dist_sums,Array4d<int> & col_dist_sums,Array4d<int> & up_col_dist_sums) const280 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>
calcDistSumsForElementInFirstRow(int i,int j,int first_col_num,Array3d<int> & dist_sums,Array4d<int> & col_dist_sums,Array4d<int> & up_col_dist_sums) const321 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