Visual Servoing Platform  version 3.6.1 under development (2024-02-13)
vpRetinex.cpp
1 /*
2  * ViSP, open source Visual Servoing Platform software.
3  * Copyright (C) 2005 - 2023 by Inria. All rights reserved.
4  *
5  * This software is free software; you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation; either version 2 of the License, or
8  * (at your option) any later version.
9  * See the file LICENSE.txt at the root directory of this source
10  * distribution for additional information about the GNU GPL.
11  *
12  * For using ViSP with software that can not be combined with the GNU
13  * GPL, please contact Inria about acquiring a ViSP Professional
14  * Edition License.
15  *
16  * See https://visp.inria.fr for more information.
17  *
18  * This software was developed at:
19  * Inria Rennes - Bretagne Atlantique
20  * Campus Universitaire de Beaulieu
21  * 35042 Rennes Cedex
22  * France
23  *
24  * If you have questions regarding the use of this file, please contact
25  * Inria at visp@inria.fr
26  *
27  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
28  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
29  *
30  * Description:
31  * Convert image types.
32  */
33 /* Retinex_.java Using ImageJ Gaussian Filter
34  * Retinex filter algorithm based on the plugin for GIMP.
35  *
36  *@author Jimenez-Hernandez Francisco <jimenezf@fi.uaemex.mx>
37  *Developed at Birmingham University, School of Dentistry. Supervised by
38  *Gabriel Landini
39  *@version 1.0
40  *
41  * 8 July 2010
42  *
43  * This version uses ImageJ Gaussian blurring instead of GIMP's linear
44  *Gaussian because there is a bug in GIMP's implementation that shifts the
45  *results of the blurring to the right of the image when using more than 3
46  *scales.
47  *
48  * Based on:
49  * MSRCR Retinex
50  * (Multi-Scale Retinex with Color Restoration)
51  * 2003 Fabien Pelisson <Fabien.Pelisson@inrialpes.fr>
52  * Retinex GIMP plug-in
53  *
54  * Copyright (C) 2009 MAO Y.B
55  * 2009. 3. 3
56  * Visual Information Processing (VIP) Group, NJUST
57  *
58  * D. J. Jobson, Z. Rahman, and G. A. Woodell. A multi-scale
59  * Retinex for bridging the gap between color images and the
60  * human observation of scenes. IEEE Transactions on Image Processing,
61  * 1997, 6(7): 965-976
62  *
63  * This program is free software; you can redistribute it and/or modify
64  * it under the terms of the GNU General Public License as published by
65  * the Free Software Foundation; either version 2 of the License, or
66  * (at your option) any later version.
67  *
68  * This program is distributed in the hope that it will be useful,
69  * but WITHOUT ANY WARRANTY; without even the implied warranty of
70  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
71  * GNU General Public License for more details.
72  *
73  * You should have received a copy of the GNU General Public License
74  * along with this program; if not, write to the Free Software
75  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
76  *
77 */
78 
84 #include <functional>
85 #include <numeric>
86 
87 #include <visp3/core/vpImageFilter.h>
88 #include <visp3/core/vpMath.h>
89 #include <visp3/imgproc/vpImgproc.h>
90 
91 #define MAX_RETINEX_SCALES 8
92 
93 std::vector<double> retinexScalesDistribution(int scaleDiv, int level, int scale)
94 {
95  std::vector<double> scales(MAX_RETINEX_SCALES);
96 
97  if (scaleDiv == 1) {
98  scales[0] = scale / 2.0;
99  }
100  else if (scaleDiv == 2) {
101  scales[0] = scale / 2.0;
102  scales[1] = scale;
103  }
104  else {
105  double size_step = scale / (double)scaleDiv;
106  int i;
107 
108  switch (level) {
109  case vp::RETINEX_UNIFORM:
110  for (i = 0; i < scaleDiv; i++) {
111  scales[(size_t)i] = 2.0 + i * size_step;
112  }
113  break;
114 
115  case vp::RETINEX_LOW:
116  size_step = std::log(scale - 2.0) / (double)scaleDiv;
117  for (i = 0; i < scaleDiv; i++) {
118  scales[(size_t)i] = 2.0 + std::pow(10.0, (i * size_step) / std::log(10.0));
119  }
120  break;
121 
122  case vp::RETINEX_HIGH:
123  size_step = std::log(scale - 2.0) / (double)scaleDiv;
124  for (i = 0; i < scaleDiv; i++) {
125  scales[(size_t)i] = scale - std::pow(10.0, (i * size_step) / std::log(10.0));
126  }
127  break;
128 
129  default:
130  break;
131  }
132  }
133 
134  return scales;
135 }
136 
137 // See: http://imagej.net/Retinex and
138 // https://docs.gimp.org/en/plug-in-retinex.html
139 void MSRCR(vpImage<vpRGBa> &I, int _scale, int scaleDiv, int level, double dynamic, int _kernelSize)
140 {
141  // Calculate the scales of filtering according to the number of filter and
142  // their distribution.
143  std::vector<double> retinexScales = retinexScalesDistribution(scaleDiv, level, _scale);
144 
145  // Filtering according to the various scales.
146  // Summarize the results of the various filters according to a specific
147  // weight(here equivalent for all).
148  double weight = 1.0 / (double)scaleDiv;
149 
150  std::vector<vpImage<double> > doubleRGB(3);
151  std::vector<vpImage<double> > doubleResRGB(3);
152  unsigned int size = I.getSize();
153 
154  int kernelSize = _kernelSize;
155  if (kernelSize == -1) {
156  // Compute the kernel size from the input image size
157  kernelSize = (int)(std::min<unsigned int>(I.getWidth(), I.getHeight()) / 2.0);
158  kernelSize = (kernelSize - kernelSize % 2) + 1;
159  }
160 
161  for (int channel = 0; channel < 3; channel++) {
162  doubleRGB[(size_t)channel] = vpImage<double>(I.getHeight(), I.getWidth());
163  doubleResRGB[(size_t)channel] = vpImage<double>(I.getHeight(), I.getWidth());
164 
165  for (unsigned int cpt = 0; cpt < size; cpt++) {
166  // Shift the pixel values by 1 to avoid problem with log(0)
167  switch (channel) {
168  case 0:
169  doubleRGB[(size_t)channel].bitmap[cpt] = I.bitmap[cpt].R + 1.0;
170  break;
171 
172  case 1:
173  doubleRGB[(size_t)channel].bitmap[cpt] = I.bitmap[cpt].G + 1.0;
174  break;
175 
176  case 2:
177  doubleRGB[(size_t)channel].bitmap[cpt] = I.bitmap[cpt].B + 1.0;
178  break;
179 
180  default:
181  break;
182  }
183  }
184 
185  for (int sc = 0; sc < scaleDiv; sc++) {
186  vpImage<double> blurImage;
187  double sigma = retinexScales[(size_t)sc];
188  vpImageFilter::gaussianBlur(doubleRGB[(size_t)channel], blurImage, (unsigned int)kernelSize, sigma);
189 
190  for (unsigned int cpt = 0; cpt < size; cpt++) {
191  // Summarize the filtered values.
192  // In fact one calculates a ratio between the original values and the
193  // filtered values.
194  doubleResRGB[(size_t)channel].bitmap[cpt] +=
195  weight * (std::log(doubleRGB[(size_t)channel].bitmap[cpt]) - std::log(blurImage.bitmap[cpt]));
196  }
197  }
198  }
199 
200  std::vector<double> dest(size * 3);
201  const double gain = 1.0, alpha = 128.0, offset = 0.0;
202 
203  for (unsigned int cpt = 0; cpt < size; cpt++) {
204  double logl = std::log((double)(I.bitmap[cpt].R + I.bitmap[cpt].G + I.bitmap[cpt].B + 3.0));
205 
206  dest[cpt * 3] = gain * (std::log(alpha * doubleRGB[0].bitmap[cpt]) - logl) * doubleResRGB[0].bitmap[cpt] + offset;
207  dest[cpt * 3 + 1] =
208  gain * (std::log(alpha * doubleRGB[1].bitmap[cpt]) - logl) * doubleResRGB[1].bitmap[cpt] + offset;
209  dest[cpt * 3 + 2] =
210  gain * (std::log(alpha * doubleRGB[2].bitmap[cpt]) - logl) * doubleResRGB[2].bitmap[cpt] + offset;
211  }
212 
213  double sum = std::accumulate(dest.begin(), dest.end(), 0.0);
214  double mean = sum / dest.size();
215 
216  std::vector<double> diff(dest.size());
217 
218 #if VISP_CXX_STANDARD > VISP_CXX_STANDARD_98
219  std::transform(dest.begin(), dest.end(), diff.begin(), std::bind(std::minus<double>(), std::placeholders::_1, mean));
220 #else
221  std::transform(dest.begin(), dest.end(), diff.begin(), std::bind2nd(std::minus<double>(), mean));
222 #endif
223 
224  double sq_sum = std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0);
225  double stdev = std::sqrt(sq_sum / dest.size());
226 
227  double mini = mean - dynamic * stdev;
228  double maxi = mean + dynamic * stdev;
229  double range = maxi - mini;
230 
231  if (vpMath::nul(range)) {
232  range = 1.0;
233  }
234 
235  for (unsigned int cpt = 0; cpt < size; cpt++) {
236  I.bitmap[cpt].R = vpMath::saturate<unsigned char>((255.0 * (dest[cpt * 3 + 0] - mini) / range));
237  I.bitmap[cpt].G = vpMath::saturate<unsigned char>((255.0 * (dest[cpt * 3 + 1] - mini) / range));
238  I.bitmap[cpt].B = vpMath::saturate<unsigned char>((255.0 * (dest[cpt * 3 + 2] - mini) / range));
239  }
240 }
241 
242 namespace vp
243 {
244 void retinex(vpImage<vpRGBa> &I, int scale, int scaleDiv, int level, const double dynamic, int kernelSize)
245 {
246  // Assert scale
247  if (scale < 16 || scale > 250) {
248  std::cerr << "Scale must be between the interval [16 - 250]" << std::endl;
249  return;
250  }
251 
252  // Assert scaleDiv
253  if (scaleDiv < 1 || scaleDiv > 8) {
254  std::cerr << "Scale division must be between the interval [1 - 8]" << std::endl;
255  return;
256  }
257 
258  if (I.getWidth() * I.getHeight() == 0) {
259  return;
260  }
261 
262  MSRCR(I, scale, scaleDiv, level, dynamic, kernelSize);
263 }
264 
265 void retinex(const vpImage<vpRGBa> &I1, vpImage<vpRGBa> &I2, int scale, int scaleDiv, int level, double dynamic,
266  int kernelSize)
267 {
268  I2 = I1;
269  vp::retinex(I2, scale, scaleDiv, level, dynamic, kernelSize);
270 }
271 };
static void gaussianBlur(const vpImage< ImageType > &I, vpImage< FilterType > &GI, unsigned int size=7, FilterType sigma=0., bool normalize=true, const vpImage< bool > *p_mask=nullptr)
unsigned int getWidth() const
Definition: vpImage.h:242
unsigned int getSize() const
Definition: vpImage.h:223
Type * bitmap
points toward the bitmap
Definition: vpImage.h:139
unsigned int getHeight() const
Definition: vpImage.h:184
static bool nul(double x, double threshold=0.001)
Definition: vpMath.h:440
VISP_EXPORT void retinex(vpImage< vpRGBa > &I, int scale=240, int scaleDiv=3, int level=RETINEX_UNIFORM, double dynamic=1.2, int kernelSize=-1)
Definition: vpRetinex.cpp:244
@ RETINEX_HIGH
Enhances the bright regions of the image.
Definition: vpImgproc.h:59
@ RETINEX_LOW
Enhances dark regions of the image.
Definition: vpImgproc.h:58
@ RETINEX_UNIFORM
Tends to treat all image intensities similarly.
Definition: vpImgproc.h:57