image_normalization.cpp
Go to the documentation of this file.
1 // *****************************************************************************
2 //
3 // Copyright (c) 2014, Southwest Research Institute® (SwRI®)
4 // All rights reserved.
5 //
6 // Redistribution and use in source and binary forms, with or without
7 // modification, are permitted provided that the following conditions are met:
8 // * Redistributions of source code must retain the above copyright
9 // notice, this list of conditions and the following disclaimer.
10 // * Redistributions in binary form must reproduce the above copyright
11 // notice, this list of conditions and the following disclaimer in the
12 // documentation and/or other materials provided with the distribution.
13 // * Neither the name of Southwest Research Institute® (SwRI®) nor the
14 // names of its contributors may be used to endorse or promote products
15 // derived from this software without specific prior written permission.
16 //
17 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 // ARE DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
21 // DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
22 // (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
23 // LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
24 // ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
26 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 //
28 // *****************************************************************************
29 /*
30  * image_normalization.cpp
31  *
32  * Created on: Jan 19, 2012
33  * Author: kkozak
34  */
35 
37 
38 namespace swri_image_util
39 {
40  void normalize_illumination(cv::Mat NormImage,
41  cv::Mat SourceImage,
42  cv::Mat& DestImage)
43  {
44  // note that Norm image should be of type CV_32FC1 and scaled between 0 and
45  // 1 (but not exactly equal to 0 anywhere)
46  cv::Mat tempIm1;
47  SourceImage.convertTo(tempIm1, CV_32FC1, 1.0, 0.0);
48  cv::divide(tempIm1, NormImage, tempIm1, 1.0);
49  cv::Mat firstImage;
50  cv::Mat secondImage;
51  tempIm1.convertTo(DestImage, CV_8UC1);
52  }
53 
54  cv::Mat generate_normalization_image(const std::vector<cv::Mat>& image_list)
55  {
56  cv::Mat norm_image;
57  if (image_list.empty())
58  {
59  return cv::Mat();
60  }
61 
62  const cv::Mat& rep_im = image_list.front();
63 
64  cv::Mat image_sum(cv::Size(rep_im.cols, rep_im.rows), CV_64F);
65 
66  for (uint32_t i = 0; i < image_list.size(); i++)
67  {
68  cv::Mat temp_im;
69  image_list[i].convertTo(temp_im, CV_64F, 1.0, 0.0);
70 
71  image_sum = image_sum + temp_im;
72  }
73 
74  image_sum = image_sum / static_cast<double>(image_list.size());
75 
76  cv::Mat mean_image;
77  image_sum.convertTo(mean_image, CV_8U);
78  cv::Mat temp_norm_image;
79 
80  cv::medianBlur(mean_image, temp_norm_image, 45);
81 
82  cv::Mat temp_norm_image2;
83  temp_norm_image.convertTo(temp_norm_image2, CV_32F);
84  double max1 = 0;
85  for (int32_t i = 0; i < temp_norm_image2.rows; i++)
86  {
87  for (int32_t j = 0; j < temp_norm_image2.cols; j++)
88  {
89  if (temp_norm_image2.at<float>(i, j) > max1)
90  {
91  max1 = temp_norm_image2.at<float>(i, j);
92  }
93  }
94  }
95 
96  temp_norm_image2 = temp_norm_image2 * (255.0 / max1);
97 
98  cv::Mat temp_norm_image3;
99  temp_norm_image2.convertTo(temp_norm_image3, CV_8U, 1.0, 0.0);
100 
101 
102  cv::GaussianBlur(temp_norm_image3,
103  norm_image,
104  cv::Size(25, 25),
105  5,
106  5);
107 
108  return norm_image;
109  }
110 
111  void MaskedBoxFilter(cv::Mat& mat, cv::Mat& mask, int32_t kernel_size)
112  {
113  mask.setTo(1, mask);
114  cv::Mat local_sums;
115  cv::boxFilter(
116  mat,
117  local_sums,
118  mat.depth(),
119  cv::Size(kernel_size, kernel_size),
120  cv::Point(-1, -1),
121  false);
122  cv::Mat local_mask_sums;
123  cv::boxFilter(
124  mask,
125  local_mask_sums,
126  mat.depth(),
127  cv::Size(kernel_size, kernel_size),
128  cv::Point(-1, -1),
129  false);
130  cv::Mat blurred;
131  cv::divide(local_sums, local_mask_sums, blurred, 1, mat.type());
132  blurred.copyTo(mat, local_mask_sums != 0);
133  }
134 
136  int32_t grid_size,
137  const cv::Mat& source_image,
138  cv::Mat& dest_image,
139  const cv::Mat& mask,
140  double max_min,
141  double min_max)
142  {
143  double cell_width = static_cast<double>(source_image.cols) / grid_size;
144  double cell_height = static_cast<double>(source_image.rows) / grid_size;
145 
146  cv::Mat max_vals(grid_size, grid_size, CV_32F);
147  cv::Mat min_vals(grid_size, grid_size, CV_32F);
148  cv::Mat grid_mask = cv::Mat::ones(grid_size, grid_size, CV_8U) * 255;
149 
150  bool has_mask = !mask.empty();
151  for(int i = 0; i < grid_size; i++)
152  {
153  for(int j = 0; j < grid_size; j++)
154  {
155  double minVal = 0;
156  double maxVal = 0;
157 
158  int x = j * cell_width;
159  int y = i * cell_height;
160  int x2 = (j + 1) * cell_width;
161  int y2 = (i + 1) * cell_height;
162  int w = x2 - x;
163  int h = y2 - y;
164 
165  cv::Rect roi = cv::Rect(x, y, w, h);
166  roi.width = std::min(source_image.cols - roi.x, roi.width);
167  roi.height = std::min(source_image.rows - roi.y, roi.height);
168 
169  if (has_mask)
170  {
171  if (cv::countNonZero(mask(roi)) > 0)
172  {
173  cv::minMaxLoc(source_image(roi), &minVal, &maxVal, 0, 0, mask(roi));
174  }
175  else
176  {
177  grid_mask.at<uint8_t>(i, j) = 0;
178  }
179  }
180  else
181  {
182  cv::minMaxLoc(source_image(roi), &minVal, &maxVal, 0, 0);
183  }
184 
185  max_vals.at<float>(i, j) = maxVal;
186  min_vals.at<float>(i, j) = minVal;
187  }
188  }
189 
190  MaskedBoxFilter(max_vals, grid_mask, 3);
191  MaskedBoxFilter(min_vals, grid_mask, 3);
192 
193  // Stretch contrast accordingly
194  for(int i = 0; i < source_image.rows; i++)
195  {
196  int cell_y = std::min(grid_size - 1, static_cast<int>(i / cell_height));
197  int cell_y_start = cell_y * cell_height;
198  int cell_y_end = std::min(source_image.rows - 1, static_cast<int>((cell_y + 1) * cell_height));
199  double py = (i - cell_y_start) / static_cast<double>(cell_y_end - cell_y_start);
200 
201  for(int j = 0; j < source_image.cols; j++)
202  {
203  // Find relevant min and max values
204  int cell_x = std::min(grid_size - 1, static_cast<int>(j / cell_width));
205  int cell_x_start = cell_x * cell_width;
206  int cell_x_end = std::min(source_image.cols - 1, static_cast<int>((cell_x + 1) * cell_width));
207  double px = (j - cell_x_start) / static_cast<double>(cell_x_end - cell_x_start);
208 
209  int cell_x1 = cell_x;
210  int cell_x2 = cell_x;
211  double px1 = 0.5;
212  double px2 = 0.5;
213  if (px < 0.5 && cell_x > 0)
214  {
215  cell_x1 = cell_x - 1;
216  cell_x2 = cell_x;
217  px2 = px + 0.5;
218  px1 = 1.0 - px2;
219  }
220  else if (px > 0.5 && cell_x + 1 < grid_size)
221  {
222  cell_x1 = cell_x;
223  cell_x2 = cell_x + 1;
224  px1 = 1.5 - px;
225  px2 = 1.0 - px1;
226  }
227 
228  int cell_y1 = cell_y;
229  int cell_y2 = cell_y;
230  double py1 = 0.5;
231  double py2 = 0.5;
232  if (py < 0.5 && cell_y > 0)
233  {
234  cell_y1 = cell_y - 1;
235  cell_y2 = cell_y;
236  py2 = py + 0.5;
237  py1 = 1.0 - py2;
238  }
239  else if (py > 0.5 && cell_y + 1 < grid_size)
240  {
241  cell_y1 = cell_y;
242  cell_y2 = cell_y + 1;
243  py1 = 1.5 - py;
244  py2 = 1.0 - py1;
245  }
246 
247  // Stretch histogram
248  double min_x1_y1 = min_vals.at<float>(cell_y1, cell_x1);
249  double max_x1_y1 = max_vals.at<float>(cell_y1, cell_x1);
250 
251  double min_x2_y1 = min_vals.at<float>(cell_y1, cell_x2);
252  double max_x2_y1 = max_vals.at<float>(cell_y1, cell_x2);
253 
254  double min_x1_y2 = min_vals.at<float>(cell_y2, cell_x1);
255  double max_x1_y2 = max_vals.at<float>(cell_y2, cell_x1);
256 
257  double min_x2_y2 = min_vals.at<float>(cell_y2, cell_x2);
258  double max_x2_y2 = max_vals.at<float>(cell_y2, cell_x2);
259 
260  //4-point interpolation
261  double xM1 = max_x1_y1 * px1 + max_x2_y1 * px2;
262  double xM2 = max_x1_y2 * px1 + max_x2_y2 * px2;
263  double M = xM1 * py1 + xM2 * py2;
264 
265  double xm1 = min_x1_y1 * px1 + min_x2_y1 * px2;
266  double xm2 = min_x1_y2 * px1 + min_x2_y2 * px2;
267  double m = xm1 * py1 + xm2 * py2;
268  double minVal = m;
269  double maxVal = M;
270 
271  if(maxVal > 255) maxVal = 255;
272  if(minVal < 0) minVal = 0;
273 
274  // Put a bound on maxVal and minVal
275  maxVal = std::max(maxVal, min_max);
276  minVal = std::min(minVal, max_min);
277 
278  double val = source_image.at<uint8_t>(i,j);
279  val = 255.0 * (val - minVal) / (maxVal - minVal);
280  if(val > 255) val = 255;
281  if(val < 0) val = 0;
282  dest_image.at<uint8_t>(i,j) = val;
283  }
284  }
285 
286  dest_image.setTo(0, mask == 0);
287  }
288 
290  const cv::Mat& src,
291  cv::Mat& dst,
292  int winsize,
293  int ftzero,
294  uchar* buf)
295  {
296  if (dst.empty())
297  {
298  dst.create(src.size(), CV_8U);
299  }
300 
301  // Source code from OpenCV: modules/calib3d/src/stereobm.cpp
302  int x, y, wsz2 = winsize / 2;
303  int* vsum = reinterpret_cast<int*>(cv::alignPtr(buf + (wsz2 + 1) * sizeof(vsum[0]), 32));
304  int scale_g = winsize * winsize / 8, scale_s = (1024 + scale_g) / (scale_g * 2);
305  const int OFS = 256 * 5, TABSZ = OFS * 2 + 256;
306  uchar tab[TABSZ];
307  const uchar* sptr = src.ptr();
308  int srcstep = src.step;
309  cv::Size size = src.size();
310 
311  scale_g *= scale_s;
312 
313  for (x = 0; x < TABSZ; x++)
314  {
315  tab[x] = (uchar)(x - OFS < -ftzero ? 0 : x - OFS > ftzero ? ftzero * 2 : x - OFS + ftzero);
316  }
317 
318  for (x = 0; x < size.width; x++)
319  {
320  vsum[x] = (ushort)(sptr[x] * (wsz2 + 2));
321  }
322 
323  for (y = 1; y < wsz2; y++)
324  {
325  for (x = 0; x < size.width; x++)
326  {
327  vsum[x] = (ushort)(vsum[x] + sptr[srcstep*y + x]);
328  }
329  }
330 
331  for (y = 0; y < size.height; y++)
332  {
333  const uchar* top = sptr + srcstep * MAX(y - wsz2 - 1, 0);
334  const uchar* bottom = sptr + srcstep * MIN(y + wsz2, size.height - 1);
335  const uchar* prev = sptr + srcstep * MAX(y - 1, 0);
336  const uchar* curr = sptr + srcstep * y;
337  const uchar* next = sptr + srcstep * MIN(y + 1, size.height - 1);
338  uchar* dptr = dst.ptr<uchar>(y);
339 
340  for (x = 0; x < size.width; x++)
341  {
342  vsum[x] = (ushort)(vsum[x] + bottom[x] - top[x]);
343  }
344 
345  for (x = 0; x <= wsz2; x++)
346  {
347  vsum[-x-1] = vsum[0];
348  vsum[size.width + x] = vsum[size.width - 1];
349  }
350 
351  int sum = vsum[0] * (wsz2 + 1);
352  for( x = 1; x <= wsz2; x++ )
353  {
354  sum += vsum[x];
355  }
356 
357  int val = ((curr[0] * 5 + curr[1] + prev[0] + next[0]) * scale_g - sum * scale_s) >> 10;
358  dptr[0] = tab[val + OFS];
359 
360  for( x = 1; x < size.width-1; x++ )
361  {
362  sum += vsum[x + wsz2] - vsum[x - wsz2 - 1];
363  val = ((curr[x] * 4 + curr[x - 1] + curr[x + 1] + prev[x] + next[x]) * scale_g - sum * scale_s) >> 10;
364  dptr[x] = tab[val + OFS];
365  }
366 
367  sum += vsum[x+wsz2] - vsum[x - wsz2 - 1];
368  val = ((curr[x] * 5 + curr[x - 1] + prev[x] + next[x]) * scale_g - sum * scale_s) >> 10;
369  dptr[x] = tab[val + OFS];
370  }
371  }
372 
373  cv::Mat scale_2_8bit(const cv::Mat& image)
374  {
375  if (image.type() == CV_8UC1)
376  return image;
377  cv::Mat Image8Bit(image.rows, image.cols, CV_8U), Image8BitColor; //Define an 8bit image
378  //Convert the image to 32bit float
379  cv::Mat ImageFloat;
380  image.convertTo(ImageFloat, CV_32F, 1, 0);
381  double maxVal; //Define the max value of image
382  cv::minMaxLoc(ImageFloat.reshape(1,1), NULL, &maxVal);//Extract the max value of image
383  //ReScale the image to 0 to 255
384  ImageFloat = ImageFloat*((1 << 8)/pow(2, ceil(log(maxVal)/log(2))));
385  ImageFloat.convertTo(Image8Bit,CV_8U, 1, 0);
386  return Image8Bit;
387  }
388 
389  cv::Mat scale_2_8bit_color(const cv::Mat& image)
390  {
391  if (image.type() == CV_8UC3)
392  return image;
393  cv::Mat Image8Bit = scale_2_8bit(image), Image8BitColor;
394  cvtColor(Image8Bit, Image8BitColor, cv::COLOR_GRAY2BGR);
395  return Image8BitColor;
396  }
397 }
TFSIMD_FORCE_INLINE const tfScalar & y() const
void NormalizeResponse(const cv::Mat &src, cv::Mat &dst, int winsize, int ftzero, uchar *buf)
cv::Mat generate_normalization_image(const std::vector< cv::Mat > &image_list)
Computes a best estimate of a normalization image from a vector of images.
CvImagePtr cvtColor(const CvImageConstPtr &source, const std::string &encoding)
void MaskedBoxFilter(cv::Mat &mat, cv::Mat &mask, int32_t kernel_size)
TFSIMD_FORCE_INLINE const tfScalar & x() const
cv::Mat scale_2_8bit_color(const cv::Mat &image)
Convert the input Mat to 8 bit color.
TFSIMD_FORCE_INLINE const tfScalar & w() const
#define NULL
void normalize_illumination(cv::Mat NormImage, cv::Mat SourceImage, cv::Mat &DestImage)
cv::Mat scale_2_8bit(const cv::Mat &image)
Convert the input Mat to 8 bit.
void ContrastStretch(int32_t grid_size, const cv::Mat &source_image, cv::Mat &dest_image, const cv::Mat &mask=cv::Mat(), double max_min=0.0, double min_max=0.0)


swri_image_util
Author(s): Kris Kozak
autogenerated on Tue Apr 6 2021 02:50:39