modules/contrib/src/retinacolor.cpp

   1 /*#******************************************************************************
   2 ** IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   3 **
   4 ** By downloading, copying, installing or using the software you agree to this license.
   5 ** If you do not agree to this license, do not download, install,
   6 ** copy or use the software.
   7 **
   8 **
   9 ** HVStools : interfaces allowing OpenCV users to integrate Human Vision System models. Presented models originate from Jeanny Herault's original research and have been reused and adapted by the author&collaborators for computed vision applications since his thesis with Alice Caplier at Gipsa-Lab.
  10 ** Use: extract still images & image sequences features, from contours details to motion spatio-temporal features, etc. for high level visual scene analysis. Also contribute to image enhancement/compression such as tone mapping.
  11 **
  12 ** Maintainers : Listic lab (code author current affiliation & applications) and Gipsa Lab (original research origins & applications)
  13 **
  14 **  Creation - enhancement process 2007-2011
  15 **      Author: Alexandre Benoit (benoit.alexandre.vision@gmail.com), LISTIC lab, Annecy le vieux, France
  16 **
  17 ** Theses algorithm have been developped by Alexandre BENOIT since his thesis with Alice Caplier at Gipsa-Lab (www.gipsa-lab.inpg.fr) and the research he pursues at LISTIC Lab (www.listic.univ-savoie.fr).
  18 ** Refer to the following research paper for more information:
  19 ** Benoit A., Caplier A., Durette B., Herault, J., "USING HUMAN VISUAL SYSTEM MODELING FOR BIO-INSPIRED LOW LEVEL IMAGE PROCESSING", Elsevier, Computer Vision and Image Understanding 114 (2010), pp. 758-773, DOI: http://dx.doi.org/10.1016/j.cviu.2010.01.011
  20 ** This work have been carried out thanks to Jeanny Herault who's research and great discussions are the basis of all this work, please take a look at his book:
  21 ** Vision: Images, Signals and Neural Networks: Models of Neural Processing in Visual Perception (Progress in Neural Processing),By: Jeanny Herault, ISBN: 9814273686. WAPI (Tower ID): 113266891.
  22 **
  23 ** The retina filter includes the research contributions of phd/research collegues from which code has been redrawn by the author :
  24 ** _take a look at the retinacolor.hpp module to discover Brice Chaix de Lavarene color mosaicing/demosaicing and the reference paper:
  25 ** ====> B. Chaix de Lavarene, D. Alleysson, B. Durette, J. Herault (2007). "Efficient demosaicing through recursive filtering", IEEE International Conference on Image Processing ICIP 2007
  26 ** _take a look at imagelogpolprojection.hpp to discover retina spatial log sampling which originates from Barthelemy Durette phd with Jeanny Herault. A Retina / V1 cortex projection is also proposed and originates from Jeanny's discussions.
  27 ** ====> more informations in the above cited Jeanny Heraults's book.
  28 **
  29 **                          License Agreement
  30 **               For Open Source Computer Vision Library
  31 **
  32 ** Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
  33 ** Copyright (C) 2008-2011, Willow Garage Inc., all rights reserved.
  34 **
  35 **               For Human Visual System tools (hvstools)
  36 ** Copyright (C) 2007-2011, LISTIC Lab, Annecy le Vieux and GIPSA Lab, Grenoble, France, all rights reserved.
  37 **
  38 ** Third party copyrights are property of their respective owners.
  39 **
  40 ** Redistribution and use in source and binary forms, with or without modification,
  41 ** are permitted provided that the following conditions are met:
  42 **
  43 ** * Redistributions of source code must retain the above copyright notice,
  44 **    this list of conditions and the following disclaimer.
  45 **
  46 ** * Redistributions in binary form must reproduce the above copyright notice,
  47 **    this list of conditions and the following disclaimer in the documentation
  48 **    and/or other materials provided with the distribution.
  49 **
  50 ** * The name of the copyright holders may not be used to endorse or promote products
  51 **    derived from this software without specific prior written permission.
  52 **
  53 ** This software is provided by the copyright holders and contributors "as is" and
  54 ** any express or implied warranties, including, but not limited to, the implied
  55 ** warranties of merchantability and fitness for a particular purpose are disclaimed.
  56 ** In no event shall the Intel Corporation or contributors be liable for any direct,
  57 ** indirect, incidental, special, exemplary, or consequential damages
  58 ** (including, but not limited to, procurement of substitute goods or services;
  59 ** loss of use, data, or profits; or business interruption) however caused
  60 ** and on any theory of liability, whether in contract, strict liability,
  61 ** or tort (including negligence or otherwise) arising in any way out of
  62 ** the use of this software, even if advised of the possibility of such damage.
  63 *******************************************************************************/
  64
  65 #include "precomp.hpp"
  66
  67 #include "retinacolor.hpp"
  68
  69 // @author Alexandre BENOIT, benoit.alexandre.vision@gmail.com, LISTIC : www.listic.univ-savoie.fr, Gipsa-Lab, France: www.gipsa-lab.inpg.fr/
  70
  71 #include <iostream>
  72 #include <ctime>
  73
  74 namespace cv
  75 {
  76
  77 // init static values
  78 static float _LMStoACr1Cr2[]={1.0,  1.0, 0.0,  1.0, -1.0, 0.0,  -0.5, -0.5, 1.0};
  79 //static double _ACr1Cr2toLMS[]={0.5,  0.5, 0.0,   0.5, -0.5, 0.0,  0.5,  0.0, 1.0};
  80 static float _LMStoLab[]={0.5774f, 0.5774f, 0.5774f, 0.4082f, 0.4082f, -0.8165f, 0.7071f, -0.7071f, 0.f};
  81
  82 // constructor/desctructor
  83 RetinaColor::RetinaColor(const unsigned int NBrows, const unsigned int NBcolumns, const RETINA_COLORSAMPLINGMETHOD samplingMethod)
  84 :BasicRetinaFilter(NBrows, NBcolumns, 3),
  85  _colorSampling(NBrows*NBcolumns),
  86  _RGBmosaic(NBrows*NBcolumns*3),
  87  _tempMultiplexedFrame(NBrows*NBcolumns),
  88  _demultiplexedTempBuffer(NBrows*NBcolumns*3),
  89  _demultiplexedColorFrame(NBrows*NBcolumns*3),
  90  _chrominance(NBrows*NBcolumns*3),
  91  _colorLocalDensity(NBrows*NBcolumns*3),
  92  _imageGradient(NBrows*NBcolumns*2)
  93 {
  94     // link to parent buffers (let's recycle !)
  95     _luminance=&_filterOutput;
  96     _multiplexedFrame=&_localBuffer;
  97
  98     _objectInit=false;
  99     _samplingMethod=samplingMethod;
 100     _saturateColors=false;
 101     _colorSaturationValue=4.0;
 102
 103     // set default spatio-temporal filter parameters
 104     setLPfilterParameters(0.0, 0.0, 1.5);
 105     setLPfilterParameters(0.0, 0.0, 10.5, 1);// for the low pass filter dedicated to contours energy extraction (demultiplexing process)
 106     setLPfilterParameters(0.f, 0.f, 0.9f, 2);
 107
 108     // init default value on image Gradient
 109     _imageGradient=0.57f;
 110
 111     // init color sampling map
 112     _initColorSampling();
 113
 114     // flush all buffers
 115     clearAllBuffers();
 116 }
 117
 118 RetinaColor::~RetinaColor()
 119 {
 120
 121 }
 122
 123 /**
 124 * function that clears all buffers of the object
 125 */
 126 void RetinaColor::clearAllBuffers()
 127 {
 128     BasicRetinaFilter::clearAllBuffers();
 129     _tempMultiplexedFrame=0.f;
 130     _demultiplexedTempBuffer=0.f;
 131
 132     _demultiplexedColorFrame=0.f;
 133     _chrominance=0.f;
 134     _imageGradient=0.57f;
 135 }
 136
 137 /**
 138 * resize retina color filter object (resize all allocated buffers)
 139 * @param NBrows: the new height size
 140 * @param NBcolumns: the new width size
 141 */
 142 void RetinaColor::resize(const unsigned int NBrows, const unsigned int NBcolumns)
 143 {
 144     BasicRetinaFilter::clearAllBuffers();
 145     _colorSampling.resize(NBrows*NBcolumns);
 146     _RGBmosaic.resize(NBrows*NBcolumns*3);
 147     _tempMultiplexedFrame.resize(NBrows*NBcolumns);
 148     _demultiplexedTempBuffer.resize(NBrows*NBcolumns*3);
 149     _demultiplexedColorFrame.resize(NBrows*NBcolumns*3);
 150     _chrominance.resize(NBrows*NBcolumns*3);
 151     _colorLocalDensity.resize(NBrows*NBcolumns*3);
 152     _imageGradient.resize(NBrows*NBcolumns*2);
 153
 154     // link to parent buffers (let's recycle !)
 155     _luminance=&_filterOutput;
 156     _multiplexedFrame=&_localBuffer;
 157
 158     // init color sampling map
 159     _initColorSampling();
 160
 161     // clean buffers
 162     clearAllBuffers();
 163 }
 164
 165
 166 void RetinaColor::_initColorSampling()
 167 {
 168
 169     // filling the conversion table for multiplexed <=> demultiplexed frame
 170     srand((unsigned)time(NULL));
 171
 172     // preInit cones probabilities
 173     _pR=_pB=_pG=0;
 174     switch (_samplingMethod)
 175     {
 176     case RETINA_COLOR_RANDOM:
 177         for (unsigned int index=0 ; index<this->getNBpixels(); ++index)
 178         {
 179
 180             // random RGB sampling
 181             unsigned int colorIndex=rand()%24;
 182
 183             if (colorIndex<8){
 184                 colorIndex=0;
 185
 186                 ++_pR;
 187             }else
 188             {
 189                 if (colorIndex<21){
 190                     colorIndex=1;
 191                     ++_pG;
 192                 }else{
 193                     colorIndex=2;
 194                     ++_pB;
 195                 }
 196             }
 197             _colorSampling[index] = colorIndex*this->getNBpixels()+index;
 198         }
 199         _pR/=(float)this->getNBpixels();
 200         _pG/=(float)this->getNBpixels();
 201         _pB/=(float)this->getNBpixels();
 202         std::cout<<"Color channels proportions: pR, pG, pB= "<<_pR<<", "<<_pG<<", "<<_pB<<", "<<std::endl;
 203         break;
 204     case RETINA_COLOR_DIAGONAL:
 205         for (unsigned int index=0 ; index<this->getNBpixels(); ++index)
 206         {
 207             _colorSampling[index] = index+((index%3+(index%_filterOutput.getNBcolumns()))%3)*_filterOutput.getNBpixels();
 208         }
 209         _pR=_pB=_pG=1.f/3;
 210         break;
 211     case RETINA_COLOR_BAYER: // default sets bayer sampling
 212         for (unsigned int index=0 ; index<_filterOutput.getNBpixels(); ++index)
 213         {
 214             //First line: R G R G
 215             _colorSampling[index] = index+((index/_filterOutput.getNBcolumns())%2)*_filterOutput.getNBpixels()+((index%_filterOutput.getNBcolumns())%2)*_filterOutput.getNBpixels();
 216             //First line: G R G R
 217             //_colorSampling[index] = 3*index+((index/_filterOutput.getNBcolumns())%2)+((index%_filterOutput.getNBcolumns()+1)%2);
 218         }
 219         _pR=_pB=0.25;
 220         _pG=0.5;
 221         break;
 222     default:
 223 #ifdef RETINACOLORDEBUG
 224         std::cerr<<"RetinaColor::No or wrong color sampling method, skeeping"<<std::endl;
 225 #endif
 226         return;
 227         break;//.. not useful, yes
 228
 229     }
 230     // feeling the mosaic buffer:
 231     _RGBmosaic=0;
 232     for (unsigned int index=0 ; index<_filterOutput.getNBpixels(); ++index)
 233         // the RGB _RGBmosaic buffer contains 1 where the pixel corresponds to a sampled color
 234         _RGBmosaic[_colorSampling[index]]=1.0;
 235
 236     // computing photoreceptors local density
 237     _spatiotemporalLPfilter(&_RGBmosaic[0], &_colorLocalDensity[0]);
 238     _spatiotemporalLPfilter(&_RGBmosaic[0]+_filterOutput.getNBpixels(), &_colorLocalDensity[0]+_filterOutput.getNBpixels());
 239     _spatiotemporalLPfilter(&_RGBmosaic[0]+_filterOutput.getDoubleNBpixels(), &_colorLocalDensity[0]+_filterOutput.getDoubleNBpixels());
 240     unsigned int maxNBpixels=3*_filterOutput.getNBpixels();
 241     register float *colorLocalDensityPTR=&_colorLocalDensity[0];
 242     for (unsigned int i=0;i<maxNBpixels;++i, ++colorLocalDensityPTR)
 243         *colorLocalDensityPTR=1.f/ *colorLocalDensityPTR;
 244
 245 #ifdef RETINACOLORDEBUG
 246     std::cout<<"INIT    _colorLocalDensity max, min: "<<_colorLocalDensity.max()<<", "<<_colorLocalDensity.min()<<std::endl;
 247 #endif
 248     // end of the init step
 249     _objectInit=true;
 250 }
 251
 252 // public functions
 253
 254 void RetinaColor::runColorDemultiplexing(const std::valarray<float> &multiplexedColorFrame, const bool adaptiveFiltering, const float maxInputValue)
 255 {
 256     // demultiplex the grey frame to RGB frame
 257     // -> first set demultiplexed frame to 0
 258     _demultiplexedTempBuffer=0;
 259     // -> demultiplex process
 260     register unsigned int *colorSamplingPRT=&_colorSampling[0];
 261     register const float *multiplexedColorFramePtr=get_data(multiplexedColorFrame);
 262     for (unsigned int indexa=0; indexa<_filterOutput.getNBpixels() ; ++indexa)
 263         _demultiplexedTempBuffer[*(colorSamplingPRT++)]=*(multiplexedColorFramePtr++);
 264
 265     // interpolate the demultiplexed frame depending on the color sampling method
 266     if (!adaptiveFiltering)
 267         _interpolateImageDemultiplexedImage(&_demultiplexedTempBuffer[0]);
 268
 269     // low pass filtering the demultiplexed frame
 270     _spatiotemporalLPfilter(&_demultiplexedTempBuffer[0], &_chrominance[0]);
 271     _spatiotemporalLPfilter(&_demultiplexedTempBuffer[0]+_filterOutput.getNBpixels(), &_chrominance[0]+_filterOutput.getNBpixels());
 272     _spatiotemporalLPfilter(&_demultiplexedTempBuffer[0]+_filterOutput.getDoubleNBpixels(), &_chrominance[0]+_filterOutput.getDoubleNBpixels());
 273
 274     /*if (_samplingMethod=BAYER)
 275     {
 276         _applyRIFfilter(_chrominance, _chrominance);
 277         _applyRIFfilter(_chrominance+_filterOutput.getNBpixels(), _chrominance+_filterOutput.getNBpixels());
 278         _applyRIFfilter(_chrominance+_filterOutput.getDoubleNBpixels(), _chrominance+_filterOutput.getDoubleNBpixels());
 279     }*/
 280
 281     // normalize by the photoreceptors local density and retrieve the local luminance
 282     register float *chrominancePTR= &_chrominance[0];
 283     register float *colorLocalDensityPTR= &_colorLocalDensity[0];
 284     register float *luminance= &(*_luminance)[0];
 285     if (!adaptiveFiltering)// compute the gradient on the luminance
 286     {
 287         if (_samplingMethod==RETINA_COLOR_RANDOM)
 288             for (unsigned int indexc=0; indexc<_filterOutput.getNBpixels() ; ++indexc, ++chrominancePTR, ++colorLocalDensityPTR, ++luminance)
 289             {
 290                 // normalize by photoreceptors density
 291                 float Cr=*(chrominancePTR)*_colorLocalDensity[indexc];
 292                 float Cg=*(chrominancePTR+_filterOutput.getNBpixels())*_colorLocalDensity[indexc+_filterOutput.getNBpixels()];
 293                 float Cb=*(chrominancePTR+_filterOutput.getDoubleNBpixels())*_colorLocalDensity[indexc+_filterOutput.getDoubleNBpixels()];
 294                 *luminance=(Cr+Cg+Cb)*_pG;
 295                 *(chrominancePTR)=Cr-*luminance;
 296                 *(chrominancePTR+_filterOutput.getNBpixels())=Cg-*luminance;
 297                 *(chrominancePTR+_filterOutput.getDoubleNBpixels())=Cb-*luminance;
 298             }
 299         else
 300             for (unsigned int indexc=0; indexc<_filterOutput.getNBpixels() ; ++indexc, ++chrominancePTR, ++colorLocalDensityPTR, ++luminance)
 301             {
 302                 float Cr=*(chrominancePTR);
 303                 float Cg=*(chrominancePTR+_filterOutput.getNBpixels());
 304                 float Cb=*(chrominancePTR+_filterOutput.getDoubleNBpixels());
 305                 *luminance=_pR*Cr+_pG*Cg+_pB*Cb;
 306                 *(chrominancePTR)=Cr-*luminance;
 307                 *(chrominancePTR+_filterOutput.getNBpixels())=Cg-*luminance;
 308                 *(chrominancePTR+_filterOutput.getDoubleNBpixels())=Cb-*luminance;
 309             }
 310
 311         // in order to get the color image, each colored map needs to be added the luminance
 312         // -> to do so, compute:  multiplexedColorFrame - remultiplexed chrominances
 313         runColorMultiplexing(_chrominance, _tempMultiplexedFrame);
 314         //lum = 1/3((f*(ImR))/(f*mR) + (f*(ImG))/(f*mG) + (f*(ImB))/(f*mB));
 315         float *luminancePTR= &(*_luminance)[0];
 316         chrominancePTR= &_chrominance[0];
 317         float *demultiplexedColorFramePTR= &_demultiplexedColorFrame[0];
 318         for (unsigned int indexp=0; indexp<_filterOutput.getNBpixels() ; ++indexp, ++luminancePTR, ++chrominancePTR, ++demultiplexedColorFramePTR)
 319         {
 320             *luminancePTR=(multiplexedColorFrame[indexp]-_tempMultiplexedFrame[indexp]);
 321             *(demultiplexedColorFramePTR)=*(chrominancePTR)+*luminancePTR;
 322             *(demultiplexedColorFramePTR+_filterOutput.getNBpixels())=*(chrominancePTR+_filterOutput.getNBpixels())+*luminancePTR;
 323             *(demultiplexedColorFramePTR+_filterOutput.getDoubleNBpixels())=*(chrominancePTR+_filterOutput.getDoubleNBpixels())+*luminancePTR;
 324         }
 325
 326     }else
 327     {
 328         register const float *multiplexedColorFramePTR= get_data(multiplexedColorFrame);
 329         for (unsigned int indexc=0; indexc<_filterOutput.getNBpixels() ; ++indexc, ++chrominancePTR, ++colorLocalDensityPTR, ++luminance, ++multiplexedColorFramePTR)
 330         {
 331             // normalize by photoreceptors density
 332             float Cr=*(chrominancePTR)*_colorLocalDensity[indexc];
 333             float Cg=*(chrominancePTR+_filterOutput.getNBpixels())*_colorLocalDensity[indexc+_filterOutput.getNBpixels()];
 334             float Cb=*(chrominancePTR+_filterOutput.getDoubleNBpixels())*_colorLocalDensity[indexc+_filterOutput.getDoubleNBpixels()];
 335             *luminance=(Cr+Cg+Cb)*_pG;
 336             _demultiplexedTempBuffer[_colorSampling[indexc]] = *multiplexedColorFramePTR - *luminance;
 337
 338         }
 339
 340         // compute the gradient of the luminance
 341         _computeGradient(&(*_luminance)[0]);
 342
 343         // adaptively filter the submosaics to get the adaptive densities, here the buffer _chrominance is used as a temp buffer
 344         _adaptiveSpatialLPfilter(&_RGBmosaic[0], &_chrominance[0]);
 345         _adaptiveSpatialLPfilter(&_RGBmosaic[0]+_filterOutput.getNBpixels(), &_chrominance[0]+_filterOutput.getNBpixels());
 346         _adaptiveSpatialLPfilter(&_RGBmosaic[0]+_filterOutput.getDoubleNBpixels(), &_chrominance[0]+_filterOutput.getDoubleNBpixels());
 347
 348         _adaptiveSpatialLPfilter(&_demultiplexedTempBuffer[0], &_demultiplexedColorFrame[0]);
 349         _adaptiveSpatialLPfilter(&_demultiplexedTempBuffer[0]+_filterOutput.getNBpixels(), &_demultiplexedColorFrame[0]+_filterOutput.getNBpixels());
 350         _adaptiveSpatialLPfilter(&_demultiplexedTempBuffer[0]+_filterOutput.getDoubleNBpixels(), &_demultiplexedColorFrame[0]+_filterOutput.getDoubleNBpixels());
 351
 352 /*      for (unsigned int index=0; index<_filterOutput.getNBpixels()*3 ; ++index) // cette boucle pourrait �tre supprimee en passant la densit� � la fonction de filtrage
 353             _demultiplexedColorFrame[index] /= _chrominance[index];*/
 354         _demultiplexedColorFrame/=_chrominance; // more optimal ;o)
 355
 356         // compute and substract the residual luminance
 357         for (unsigned int index=0; index<_filterOutput.getNBpixels() ; ++index)
 358         {
 359             float residu = _pR*_demultiplexedColorFrame[index] + _pG*_demultiplexedColorFrame[index+_filterOutput.getNBpixels()] + _pB*_demultiplexedColorFrame[index+_filterOutput.getDoubleNBpixels()];
 360             _demultiplexedColorFrame[index] = _demultiplexedColorFrame[index] - residu;
 361             _demultiplexedColorFrame[index+_filterOutput.getNBpixels()] = _demultiplexedColorFrame[index+_filterOutput.getNBpixels()] - residu;
 362             _demultiplexedColorFrame[index+_filterOutput.getDoubleNBpixels()] = _demultiplexedColorFrame[index+_filterOutput.getDoubleNBpixels()] - residu;
 363         }
 364
 365         // multiplex the obtained chrominance
 366         runColorMultiplexing(_demultiplexedColorFrame, _tempMultiplexedFrame);
 367         _demultiplexedTempBuffer=0;
 368
 369         // get the luminance, et and add it to each chrominance
 370         for (unsigned int index=0; index<_filterOutput.getNBpixels() ; ++index)
 371         {
 372             (*_luminance)[index]=multiplexedColorFrame[index]-_tempMultiplexedFrame[index];
 373             _demultiplexedTempBuffer[_colorSampling[index]] = _demultiplexedColorFrame[_colorSampling[index]];//multiplexedColorFrame[index] - (*_luminance)[index];
 374         }
 375
 376         _spatiotemporalLPfilter(&_demultiplexedTempBuffer[0], &_demultiplexedTempBuffer[0]);
 377         _spatiotemporalLPfilter(&_demultiplexedTempBuffer[0]+_filterOutput.getNBpixels(), &_demultiplexedTempBuffer[0]+_filterOutput.getNBpixels());
 378         _spatiotemporalLPfilter(&_demultiplexedTempBuffer[0]+_filterOutput.getDoubleNBpixels(), &_demultiplexedTempBuffer[0]+_filterOutput.getDoubleNBpixels());
 379
 380         // get the luminance and add it to each chrominance
 381         for (unsigned int index=0; index<_filterOutput.getNBpixels() ; ++index)
 382         {
 383             _demultiplexedColorFrame[index] = _demultiplexedTempBuffer[index]*_colorLocalDensity[index]+ (*_luminance)[index];
 384             _demultiplexedColorFrame[index+_filterOutput.getNBpixels()] = _demultiplexedTempBuffer[index+_filterOutput.getNBpixels()]*_colorLocalDensity[index+_filterOutput.getNBpixels()]+ (*_luminance)[index];
 385             _demultiplexedColorFrame[index+_filterOutput.getDoubleNBpixels()] = _demultiplexedTempBuffer[index+_filterOutput.getDoubleNBpixels()]*_colorLocalDensity[index+_filterOutput.getDoubleNBpixels()]+ (*_luminance)[index];
 386         }
 387     }
 388
 389     // eliminate saturated colors by simple clipping values to the input range
 390     clipRGBOutput_0_maxInputValue(NULL, maxInputValue);
 391
 392     /* transfert image gradient in order to check validity
 393     memcpy((*_luminance), _imageGradient, sizeof(float)*_filterOutput.getNBpixels());
 394     memcpy(_demultiplexedColorFrame, _imageGradient+_filterOutput.getNBpixels(), sizeof(float)*_filterOutput.getNBpixels());
 395     memcpy(_demultiplexedColorFrame+_filterOutput.getNBpixels(), _imageGradient+_filterOutput.getNBpixels(), sizeof(float)*_filterOutput.getNBpixels());
 396     memcpy(_demultiplexedColorFrame+2*_filterOutput.getNBpixels(), _imageGradient+_filterOutput.getNBpixels(), sizeof(float)*_filterOutput.getNBpixels());
 397      */
 398
 399     if (_saturateColors)
 400     {
 401         TemplateBuffer<float>::normalizeGrayOutputCentredSigmoide(128, _colorSaturationValue, maxInputValue, &_demultiplexedColorFrame[0], &_demultiplexedColorFrame[0], _filterOutput.getNBpixels());
 402         TemplateBuffer<float>::normalizeGrayOutputCentredSigmoide(128, _colorSaturationValue, maxInputValue, &_demultiplexedColorFrame[0]+_filterOutput.getNBpixels(), &_demultiplexedColorFrame[0]+_filterOutput.getNBpixels(), _filterOutput.getNBpixels());
 403         TemplateBuffer<float>::normalizeGrayOutputCentredSigmoide(128, _colorSaturationValue, maxInputValue, &_demultiplexedColorFrame[0]+_filterOutput.getNBpixels()*2, &_demultiplexedColorFrame[0]+_filterOutput.getNBpixels()*2, _filterOutput.getNBpixels());
 404     }
 405 }
 406
 407 // color multiplexing: input frame size=_NBrows*_filterOutput.getNBcolumns()*3, multiplexedFrame output size=_NBrows*_filterOutput.getNBcolumns()
 408 void RetinaColor::runColorMultiplexing(const std::valarray<float> &demultiplexedInputFrame, std::valarray<float> &multiplexedFrame)
 409 {
 410     // multiply each color layer by its bayer mask
 411     register unsigned int *colorSamplingPTR= &_colorSampling[0];
 412     register float *multiplexedFramePTR= &multiplexedFrame[0];
 413     for (unsigned int indexp=0; indexp<_filterOutput.getNBpixels(); ++indexp)
 414         *(multiplexedFramePTR++)=demultiplexedInputFrame[*(colorSamplingPTR++)];
 415 }
 416
 417 void RetinaColor::normalizeRGBOutput_0_maxOutputValue(const float maxOutputValue)
 418 {
 419     //normalizeGrayOutputCentredSigmoide(0.0, 2, _chrominance);
 420     TemplateBuffer<float>::normalizeGrayOutput_0_maxOutputValue(&_demultiplexedColorFrame[0], 3*_filterOutput.getNBpixels(), maxOutputValue);
 421     //normalizeGrayOutputCentredSigmoide(0.0, 2, _chrominance+_filterOutput.getNBpixels());
 422     //normalizeGrayOutput_0_maxOutputValue(_demultiplexedColorFrame+_filterOutput.getNBpixels(), _filterOutput.getNBpixels(), maxOutputValue);
 423     //normalizeGrayOutputCentredSigmoide(0.0, 2, _chrominance+2*_filterOutput.getNBpixels());
 424     //normalizeGrayOutput_0_maxOutputValue(_demultiplexedColorFrame+_filterOutput.getDoubleNBpixels(), _filterOutput.getNBpixels(), maxOutputValue);
 425     TemplateBuffer<float>::normalizeGrayOutput_0_maxOutputValue(&(*_luminance)[0], _filterOutput.getNBpixels(), maxOutputValue);
 426 }
 427
 428 /// normalize output between 0 and maxOutputValue;
 429 void RetinaColor::clipRGBOutput_0_maxInputValue(float *inputOutputBuffer, const float maxInputValue)
 430 {
 431     //std::cout<<"RetinaColor::normalizing RGB frame..."<<std::endl;
 432     // if outputBuffer unsassigned, the rewrite the buffer
 433     if (inputOutputBuffer==NULL)
 434         inputOutputBuffer= &_demultiplexedColorFrame[0];
 435
 436 #ifdef MAKE_PARALLEL // call the TemplateBuffer TBB clipping method
 437         cv::parallel_for_(cv::Range(0,_filterOutput.getNBpixels()*3), Parallel_clipBufferValues<float>(inputOutputBuffer, 0,  maxInputValue));
 438 #else
 439     register float *inputOutputBufferPTR=inputOutputBuffer;
 440     for (register unsigned int jf = 0; jf < _filterOutput.getNBpixels()*3; ++jf, ++inputOutputBufferPTR)
 441     {
 442         if (*inputOutputBufferPTR>maxInputValue)
 443             *inputOutputBufferPTR=maxInputValue;
 444         else if (*inputOutputBufferPTR<0)
 445             *inputOutputBufferPTR=0;
 446     }
 447 #endif
 448     //std::cout<<"RetinaColor::...normalizing RGB frame OK"<<std::endl;
 449 }
 450
 451 void RetinaColor::_interpolateImageDemultiplexedImage(float *inputOutputBuffer)
 452 {
 453
 454     switch(_samplingMethod)
 455     {
 456
 457     case RETINA_COLOR_RANDOM:
 458         return; // no need to interpolate
 459         break;
 460
 461     case RETINA_COLOR_DIAGONAL:
 462         _interpolateSingleChannelImage111(inputOutputBuffer);
 463         break;
 464
 465     case RETINA_COLOR_BAYER: // default sets bayer sampling
 466         _interpolateBayerRGBchannels(inputOutputBuffer);
 467         break;
 468     default:
 469         std::cerr<<"RetinaColor::No or wrong color sampling method, skeeping"<<std::endl;
 470         return;
 471         break;//.. not useful, yes
 472
 473     }
 474
 475 }
 476
 477 void RetinaColor::_interpolateSingleChannelImage111(float *inputOutputBuffer)
 478 {
 479     for (unsigned int indexr=0 ; indexr<_filterOutput.getNBrows(); ++indexr)
 480     {
 481         for (unsigned int indexc=1 ; indexc<_filterOutput.getNBcolumns()-1; ++indexc)
 482         {
 483             unsigned int index=indexc+indexr*_filterOutput.getNBcolumns();
 484             inputOutputBuffer[index]=(inputOutputBuffer[index-1]+inputOutputBuffer[index]+inputOutputBuffer[index+1])/3.f;
 485         }
 486     }
 487     for (unsigned int indexc=0 ; indexc<_filterOutput.getNBcolumns(); ++indexc)
 488     {
 489         for (unsigned int indexr=1 ; indexr<_filterOutput.getNBrows()-1; ++indexr)
 490         {
 491             unsigned int index=indexc+indexr*_filterOutput.getNBcolumns();
 492             inputOutputBuffer[index]=(inputOutputBuffer[index-_filterOutput.getNBcolumns()]+inputOutputBuffer[index]+inputOutputBuffer[index+_filterOutput.getNBcolumns()])/3.f;
 493         }
 494     }
 495 }
 496
 497 void RetinaColor::_interpolateBayerRGBchannels(float *inputOutputBuffer)
 498 {
 499     for (unsigned int indexr=0 ; indexr<_filterOutput.getNBrows()-1; indexr+=2)
 500     {
 501         for (unsigned int indexc=1 ; indexc<_filterOutput.getNBcolumns()-1; indexc+=2)
 502         {
 503             unsigned int indexR=indexc+indexr*_filterOutput.getNBcolumns();
 504             unsigned int indexB=_filterOutput.getDoubleNBpixels()+indexc+1+(indexr+1)*_filterOutput.getNBcolumns();
 505             inputOutputBuffer[indexR]=(inputOutputBuffer[indexR-1]+inputOutputBuffer[indexR+1])/2.f;
 506             inputOutputBuffer[indexB]=(inputOutputBuffer[indexB-1]+inputOutputBuffer[indexB+1])/2.f;
 507         }
 508     }
 509     for (unsigned int indexr=1 ; indexr<_filterOutput.getNBrows()-1; indexr+=2)
 510     {
 511         for (unsigned int indexc=0 ; indexc<_filterOutput.getNBcolumns(); ++indexc)
 512         {
 513             unsigned int indexR=indexc+indexr*_filterOutput.getNBcolumns();
 514             unsigned int indexB=_filterOutput.getDoubleNBpixels()+indexc+1+(indexr+1)*_filterOutput.getNBcolumns();
 515             inputOutputBuffer[indexR]=(inputOutputBuffer[indexR-_filterOutput.getNBcolumns()]+inputOutputBuffer[indexR+_filterOutput.getNBcolumns()])/2.f;
 516             inputOutputBuffer[indexB]=(inputOutputBuffer[indexB-_filterOutput.getNBcolumns()]+inputOutputBuffer[indexB+_filterOutput.getNBcolumns()])/2.f;
 517
 518         }
 519     }
 520     for (unsigned int indexr=1 ; indexr<_filterOutput.getNBrows()-1; ++indexr)
 521         for (unsigned int indexc=0 ; indexc<_filterOutput.getNBcolumns(); indexc+=2)
 522         {
 523             unsigned int indexG=_filterOutput.getNBpixels()+indexc+(indexr)*_filterOutput.getNBcolumns()+indexr%2;
 524             inputOutputBuffer[indexG]=(inputOutputBuffer[indexG-1]+inputOutputBuffer[indexG+1]+inputOutputBuffer[indexG-_filterOutput.getNBcolumns()]+inputOutputBuffer[indexG+_filterOutput.getNBcolumns()])*0.25f;
 525         }
 526 }
 527
 528 void RetinaColor::_applyRIFfilter(const float *sourceBuffer, float *destinationBuffer)
 529 {
 530     for (unsigned int indexr=1 ; indexr<_filterOutput.getNBrows()-1; ++indexr)
 531     {
 532         for (unsigned int indexc=1 ; indexc<_filterOutput.getNBcolumns()-1; ++indexc)
 533         {
 534             unsigned int index=indexc+indexr*_filterOutput.getNBcolumns();
 535             _tempMultiplexedFrame[index]=(4.f*sourceBuffer[index]+sourceBuffer[index-1-_filterOutput.getNBcolumns()]+sourceBuffer[index-1+_filterOutput.getNBcolumns()]+sourceBuffer[index+1-_filterOutput.getNBcolumns()]+sourceBuffer[index+1+_filterOutput.getNBcolumns()])*0.125f;
 536         }
 537     }
 538     memcpy(destinationBuffer, &_tempMultiplexedFrame[0], sizeof(float)*_filterOutput.getNBpixels());
 539 }
 540
 541 void RetinaColor::_getNormalizedContoursImage(const float *inputFrame, float *outputFrame)
 542 {
 543     float maxValue=0.f;
 544     float normalisationFactor=1.f/3.f;
 545     for (unsigned int indexr=1 ; indexr<_filterOutput.getNBrows()-1; ++indexr)
 546     {
 547         for (unsigned int indexc=1 ; indexc<_filterOutput.getNBcolumns()-1; ++indexc)
 548         {
 549             unsigned int index=indexc+indexr*_filterOutput.getNBcolumns();
 550             outputFrame[index]=normalisationFactor*fabs(8.f*inputFrame[index]-inputFrame[index-1]-inputFrame[index+1]-inputFrame[index-_filterOutput.getNBcolumns()]-inputFrame[index+_filterOutput.getNBcolumns()]-inputFrame[index-1-_filterOutput.getNBcolumns()]-inputFrame[index-1+_filterOutput.getNBcolumns()]-inputFrame[index+1-_filterOutput.getNBcolumns()]-inputFrame[index+1+_filterOutput.getNBcolumns()]);
 551             if (outputFrame[index]>maxValue)
 552                 maxValue=outputFrame[index];
 553         }
 554     }
 555     normalisationFactor=1.f/maxValue;
 556     // normalisation [0, 1]
 557     for (unsigned int indexp=1 ; indexp<_filterOutput.getNBrows()-1; ++indexp)
 558        outputFrame[indexp]=outputFrame[indexp]*normalisationFactor;
 559 }
 560
 561 //////////////////////////////////////////////////////////
 562 //        ADAPTIVE BASIC RETINA FILTER
 563 //////////////////////////////////////////////////////////
 564 // run LP filter for a new frame input and save result at a specific output adress
 565 void RetinaColor::_adaptiveSpatialLPfilter(const float *inputFrame, float *outputFrame)
 566 {
 567
 568     /**********/
 569     _gain = (1-0.57f)*(1-0.57f)*(1-0.06f)*(1-0.06f);
 570
 571     // launch the serie of 1D directional filters in order to compute the 2D low pass filter
 572     // -> horizontal filters work with the first layer of imageGradient
 573     _adaptiveHorizontalCausalFilter_addInput(inputFrame, outputFrame, 0, _filterOutput.getNBrows());
 574     _horizontalAnticausalFilter_Irregular(outputFrame, 0, _filterOutput.getNBrows(), &_imageGradient[0]);
 575     // -> horizontal filters work with the second layer of imageGradient
 576     _verticalCausalFilter_Irregular(outputFrame, 0, _filterOutput.getNBcolumns(), &_imageGradient[0]+_filterOutput.getNBpixels());
 577     _adaptiveVerticalAnticausalFilter_multGain(outputFrame, 0, _filterOutput.getNBcolumns());
 578 }
 579
 580 //  horizontal causal filter which adds the input inside... replaces the parent _horizontalCausalFilter_Irregular_addInput by avoiding a product for each pixel
 581 void RetinaColor::_adaptiveHorizontalCausalFilter_addInput(const float *inputFrame, float *outputFrame, unsigned int IDrowStart, unsigned int IDrowEnd)
 582 {
 583 #ifdef MAKE_PARALLEL
 584         cv::parallel_for_(cv::Range(IDrowStart,IDrowEnd), Parallel_adaptiveHorizontalCausalFilter_addInput(inputFrame, outputFrame, &_imageGradient[0], _filterOutput.getNBcolumns()));
 585 #else
 586     register float* outputPTR=outputFrame+IDrowStart*_filterOutput.getNBcolumns();
 587     register const float* inputPTR=inputFrame+IDrowStart*_filterOutput.getNBcolumns();
 588     register const float *imageGradientPTR= &_imageGradient[0]+IDrowStart*_filterOutput.getNBcolumns();
 589     for (unsigned int IDrow=IDrowStart; IDrow<IDrowEnd; ++IDrow)
 590     {
 591         register float result=0;
 592         for (unsigned int index=0; index<_filterOutput.getNBcolumns(); ++index)
 593         {
 594             //std::cout<<(*imageGradientPTR)<<" ";
 595             result = *(inputPTR++) + (*imageGradientPTR)* result;
 596             *(outputPTR++) = result;
 597             ++imageGradientPTR;
 598         }
 599         //        std::cout<<" "<<std::endl;
 600     }
 601 #endif
 602 }
 603
 604 //  vertical anticausal filter which multiplies the output by _gain... replaces the parent _verticalAnticausalFilter_multGain by avoiding a product for each pixel and taking into account the second layer of the _imageGradient buffer
 605 void RetinaColor::_adaptiveVerticalAnticausalFilter_multGain(float *outputFrame, unsigned int IDcolumnStart, unsigned int IDcolumnEnd)
 606 {
 607 #ifdef MAKE_PARALLEL
 608         cv::parallel_for_(cv::Range(IDcolumnStart,IDcolumnEnd), Parallel_adaptiveVerticalAnticausalFilter_multGain(outputFrame, &_imageGradient[0]+_filterOutput.getNBpixels(), _filterOutput.getNBrows(), _filterOutput.getNBcolumns(), _gain));
 609 #else
 610     float* outputOffset=outputFrame+_filterOutput.getNBpixels()-_filterOutput.getNBcolumns();
 611     float* gradOffset= &_imageGradient[0]+_filterOutput.getNBpixels()*2-_filterOutput.getNBcolumns();
 612
 613     for (unsigned int IDcolumn=IDcolumnStart; IDcolumn<IDcolumnEnd; ++IDcolumn)
 614     {
 615         register float result=0;
 616         register float *outputPTR=outputOffset+IDcolumn;
 617         register float *imageGradientPTR=gradOffset+IDcolumn;
 618         for (unsigned int index=0; index<_filterOutput.getNBrows(); ++index)
 619         {
 620             result = *(outputPTR) + (*(imageGradientPTR)) * result;
 621             *(outputPTR) = _gain*result;
 622             outputPTR-=_filterOutput.getNBcolumns();
 623             imageGradientPTR-=_filterOutput.getNBcolumns();
 624         }
 625     }
 626 #endif
 627 }
 628
 629 ///////////////////////////
 630 void RetinaColor::_computeGradient(const float *luminance)
 631 {
 632     for (unsigned int idLine=2;idLine<_filterOutput.getNBrows()-2;++idLine)
 633     {
 634         for (unsigned int idColumn=2;idColumn<_filterOutput.getNBcolumns()-2;++idColumn)
 635         {
 636             const unsigned int pixelIndex=idColumn+_filterOutput.getNBcolumns()*idLine;
 637
 638             // horizontal and vertical local gradients
 639             const float verticalGrad=fabs(luminance[pixelIndex+_filterOutput.getNBcolumns()]-luminance[pixelIndex-_filterOutput.getNBcolumns()]);
 640             const float horizontalGrad=fabs(luminance[pixelIndex+1]-luminance[pixelIndex-1]);
 641
 642             // neighborhood horizontal and vertical gradients
 643             const float verticalGrad_p=fabs(luminance[pixelIndex]-luminance[pixelIndex-2*_filterOutput.getNBcolumns()]);
 644             const float horizontalGrad_p=fabs(luminance[pixelIndex]-luminance[pixelIndex-2]);
 645             const float verticalGrad_n=fabs(luminance[pixelIndex+2*_filterOutput.getNBcolumns()]-luminance[pixelIndex]);
 646             const float horizontalGrad_n=fabs(luminance[pixelIndex+2]-luminance[pixelIndex]);
 647
 648             const float horizontalGradient=0.5f*horizontalGrad+0.25f*(horizontalGrad_p+horizontalGrad_n);
 649             const float verticalGradient=0.5f*verticalGrad+0.25f*(verticalGrad_p+verticalGrad_n);
 650
 651             // compare local gradient means and fill the appropriate filtering coefficient value that will be used in adaptative filters
 652             if (horizontalGradient<verticalGradient)
 653             {
 654                 _imageGradient[pixelIndex+_filterOutput.getNBpixels()]=0.06f;
 655                 _imageGradient[pixelIndex]=0.57f;
 656             }
 657             else
 658             {
 659                 _imageGradient[pixelIndex+_filterOutput.getNBpixels()]=0.57f;
 660                 _imageGradient[pixelIndex]=0.06f;
 661             }
 662         }
 663     }
 664 }
 665
 666 bool RetinaColor::applyKrauskopfLMS2Acr1cr2Transform(std::valarray<float> &result)
 667 {
 668     bool processSuccess=true;
 669     // basic preliminary error check
 670     if (result.size()!=_demultiplexedColorFrame.size())
 671     {
 672         std::cerr<<"RetinaColor::applyKrauskopfLMS2Acr1cr2Transform: input buffer does not match retina buffer size, conversion aborted"<<std::endl;
 673         return false;
 674     }
 675
 676     // apply transformation
 677     _applyImageColorSpaceConversion(_demultiplexedColorFrame, result, _LMStoACr1Cr2);
 678
 679     return processSuccess;
 680 }
 681
 682 bool RetinaColor::applyLMS2LabTransform(std::valarray<float> &result)
 683 {
 684     bool processSuccess=true;
 685     // basic preliminary error check
 686     if (result.size()!=_demultiplexedColorFrame.size())
 687     {
 688         std::cerr<<"RetinaColor::applyKrauskopfLMS2Acr1cr2Transform: input buffer does not match retina buffer size, conversion aborted"<<std::endl;
 689         return false;
 690     }
 691
 692     // apply transformation
 693     _applyImageColorSpaceConversion(_demultiplexedColorFrame, result, _LMStoLab);
 694
 695     return processSuccess;
 696 }
 697
 698 // template function able to perform a custom color space transformation
 699 void RetinaColor::_applyImageColorSpaceConversion(const std::valarray<float> &inputFrameBuffer, std::valarray<float> &outputFrameBuffer, const float *transformTable)
 700 {
 701     // two step methods in order to allow inputFrame and outputFrame to be the same
 702     unsigned int nbPixels=(unsigned int)(inputFrameBuffer.size()/3), dbpixels=(unsigned int)(2*inputFrameBuffer.size()/3);
 703
 704     const float *inputFrame=get_data(inputFrameBuffer);
 705     float *outputFrame= &outputFrameBuffer[0];
 706
 707     for (unsigned int dataIndex=0; dataIndex<nbPixels;++dataIndex, ++outputFrame, ++inputFrame)
 708     {
 709         // first step, compute each new values
 710         float layer1 = *(inputFrame)**(transformTable+0)  +*(inputFrame+nbPixels)**(transformTable+1)  +*(inputFrame+dbpixels)**(transformTable+2);
 711         float layer2 = *(inputFrame)**(transformTable+3)  +*(inputFrame+nbPixels)**(transformTable+4)  +*(inputFrame+dbpixels)**(transformTable+5);
 712         float layer3 = *(inputFrame)**(transformTable+6)  +*(inputFrame+nbPixels)**(transformTable+7)  +*(inputFrame+dbpixels)**(transformTable+8);
 713         // second, affect the output
 714         *(outputFrame)          = layer1;
 715         *(outputFrame+nbPixels) = layer2;
 716         *(outputFrame+dbpixels) = layer3;
 717     }
 718 }
 719
 720 }