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