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
getNearestPowerOf2(int value)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>
FastNlMeansDenoisingInvoker(const Mat & src,Mat & dst,int template_window_size,int search_window_size,const float * h)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>
operator ()(const Range & range) const149 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>
calcDistSumsForFirstElementInRow(int i,Array2d<int> & dist_sums,Array3d<int> & col_dist_sums,Array3d<int> & up_col_dist_sums) const253 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>
calcDistSumsForElementInFirstRow(int i,int j,int first_col_num,Array2d<int> & dist_sums,Array3d<int> & col_dist_sums,Array3d<int> & up_col_dist_sums) const287 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