Visual Servoing Platform version 3.5.0
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
98std::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) {
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
142void 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
264void 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
306void 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}
static void gaussianBlur(const vpImage< unsigned char > &I, vpImage< double > &GI, unsigned int size=7, double sigma=0., bool normalize=true)
unsigned int getWidth() const
Definition: vpImage.h:246
unsigned int getSize() const
Definition: vpImage.h:227
Type * bitmap
points toward the bitmap
Definition: vpImage.h:143
unsigned int getHeight() const
Definition: vpImage.h:188
static bool nul(double x, double s=0.001)
Definition: vpMath.h:286
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
@ RETINEX_HIGH
Definition: vpImgproc.h:56
@ RETINEX_LOW
Definition: vpImgproc.h:56
@ RETINEX_UNIFORM
Definition: vpImgproc.h:56