Visual Servoing Platform version 3.6.0
Loading...
Searching...
No Matches
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
93std::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) {
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
139void 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(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
242namespace vp
243{
244void 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
265void 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< unsigned char > &I, vpImage< FilterType > &GI, unsigned int size=7, FilterType sigma=0., bool normalize=true)
Definition of the vpImage class member functions.
Definition vpImage.h:135
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:360
static _Tp saturate(unsigned char v)
Definition vpMath.h:221
VISP_EXPORT void retinex(vpImage< vpRGBa > &I, int scale=240, int scaleDiv=3, int level=RETINEX_UNIFORM, double dynamic=1.2, int kernelSize=-1)
@ 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