HDRip/HDRip/ImageHDR.cpp

// This file is part of HDRip.
// 
//     HDRip is free software: you can redistribute it and/or modify it
//     under the terms of the GNU General Public License as published by
//     the Free Software Foundation, either version 3 of the License, or
//     (at your option) any later version.
// 
//     HDRip is distributed in the hope that it will be useful, but
//     WITHOUT ANY WARRANTY; without even the implied warranty of
//     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//     GNU General Public License for more details.
// 
//     You should have received a copy of the GNU General Public License
//     along with HDRip.   If not, see <https://www.gnu.org/licenses/>.
//
// HDRip project
// Author : Rémi Synave
// Contact : remi.synave@univ-littoral.fr

#include "pch.h"

#include "ImageHDR.hpp"
#include "MT_exposure.hpp"
#include "MT_contrast.hpp"
#include "MT_histogram_regularization.hpp"
#include "MT_lightnessMask.hpp"
#include "MT_saturation.hpp"
#include "MT_colorEditor.hpp"
#include "Conversion.hpp"
#include "YCurve.hpp"
#include "Utils.hpp"

#include <cmath>
#include <thread>
#include <iostream>
#include <vector>
#include <ctime>
#include <cstdlib>


/* Member methods*/


ImageHDR::ImageHDR(float* d, unsigned int w, unsigned int h)
{
  width = w;
  height = h;
  this->min_intensity = d[0];
  this->max_intensity = d[0];
  data = new float[width*height*3];
  for (unsigned int i = 0; i < width*height*3; i++) {
    data[i]=d[i];
    if (this->min_intensity > data[i])
      this->min_intensity = data[i];
    if (this->max_intensity < data[i])
      this->max_intensity = data[i];
  }
  linear = true;
  colorspace = Colorspace::RGB;
}


void ImageHDR::display_pixel(unsigned int i) const
{
  if (linear)
    std::cout << "LINEAIRE - ";
  else
    std::cout << "NON LINEAIRE - ";

  std::cout << "Pixel : ( " << data[i * 3] << "   " << data[i * 3 + 1] << "   " << data[i * 3 + 2] << " )" << std::endl;
}


void ImageHDR::display_pixel(unsigned int i, unsigned int j) const
{
  display_pixel(j * width + i);
}



/****************************************/
/**************** LINEAR ****************/
/****************************************/

void ImageHDR::linear_to_non_linear()
{
  float* non_linear = Conversion::linear_to_non_linear(data, width * height * 3);
  delete[](data);
  data = non_linear;
}


void ImageHDR::non_linear_to_linear()
{
  float* linear = Conversion::non_linear_to_linear(data, width * height * 3);
  delete[](data);
  data = linear;
}



/****************************************/
/*************** EXPOSURE ***************/
/****************************************/

#ifdef _MT_

void* exposure_MT(void* arg)
{
  MT_exposure* a = (MT_exposure*)arg;

  float* data = a->data;

  for (unsigned int i = 0; i < a->length; i++)
    data[i] *= a->coeff;

  return arg;
}

void ImageHDR::exposure(const float ev)
{
  float coeff = powf(2, ev);

  if (!linear)
    {
      non_linear_to_linear();
      linear = true;
    }

  std::thread tab_t[_MT_];
  MT_exposure tab_a[_MT_];

  unsigned int id;
  unsigned int tab_length = width * height * 3;
  unsigned int block_size = tab_length / _MT_;
  for (id = 0; id < _MT_; id++) {
    tab_a[id].data = data + (id * block_size);
    tab_a[id].length = block_size;
    tab_a[id].coeff = coeff;
    
    if (id == (_MT_ - 1))
      tab_a[id].length = tab_length - ((_MT_ - 1) * block_size);

    tab_t[id] = std::thread(exposure_MT, (void*)(tab_a + id));
  }

  for (id = 0; id < _MT_; id++) {
    tab_t[id].join();
  }

}

#else

void ImageHDR::exposure(const float ev)
{
  float coef = powf(2, ev);

  if (!linear)
    {
      non_linear_to_linear();
      linear = true;
    }

  for (unsigned int i = 0; i < width * height * 3; i++)
    data[i] *= coef;
}

#endif



/****************************************/
/*************** CONTRAST ***************/
/****************************************/

#ifdef _MT_

void* contrast_MT(void* arg)
{
  MT_contrast* a = (MT_contrast*)arg;

  float* data = a->data;

  for (unsigned int i = 0; i < a->length; i++)
    data[i] = a->coeff * (data[i] - 0.5f) + 0.5f;

  return arg;
}

void ImageHDR::contrast(const float c)
{
  float max_contrast_factor = 2.0f, scaling_factor = 1.0f, contrast_value = c;

  if (linear)
    {
      linear_to_non_linear();
      linear = false;
    }

  if (contrast_value != 0.0f)
    {
      contrast_value = contrast_value / 100.0f;
      if (contrast_value > 0.0f)
	{
	  scaling_factor = 1 * (1 - contrast_value) + max_contrast_factor * contrast_value;
	}
      else
	{
	  contrast_value = -contrast_value;
	  scaling_factor = 1 * (1 - contrast_value) + max_contrast_factor * contrast_value;
	  scaling_factor = 1 / scaling_factor;
	}
    }

  std::thread tab_t[_MT_];
  MT_contrast tab_a[_MT_];

  unsigned int id;
  unsigned int tab_length = width * height * 3;
  unsigned int block_size = tab_length / _MT_;
  for (id = 0; id < _MT_; id++)
    {
      tab_a[id].data = data + (id * block_size);
      tab_a[id].length = block_size;
      tab_a[id].coeff = scaling_factor;

      if (id == (_MT_ - 1))
	tab_a[id].length = tab_length - ((_MT_ - 1) * block_size);

      tab_t[id] = std::thread(contrast_MT, (void*)(tab_a + id));
    }

  for (id = 0; id < _MT_; id++) {
    tab_t[id].join();
  }
}

#else

void ImageHDR::contrast(const float c)
{
  float max_contrast_factor = 2.0f, scaling_factor = 1.0f, contrast_value = c;

  if (linear)
    {
      linear_to_non_linear();
      linear = false;
    }

  if (contrast_value != 0.0f)
    {
      contrast_value = contrast_value / 100.0f;
      if (contrast_value > 0.0f)
	{
	  scaling_factor = 1 * (1 - contrast_value) + max_contrast_factor * contrast_value;
	}
      else
	{
	  contrast_value = -contrast_value;
	  scaling_factor = 1 * (1 - contrast_value) + max_contrast_factor * contrast_value;
	  scaling_factor = 1 / scaling_factor;
	}
    }

  for (unsigned int i = 0; i < width * height * 3; i++)
    data[i] = scaling_factor * (data[i] - 0.5f) + 0.5f;

}

#endif



/****************************************/
/**************** YCURVE ****************/
/****************************************/

#ifdef _MT_

void* histogram_regularization_MT(void* arg)
{
  MT_histogram_regularization* a = (MT_histogram_regularization*)arg;

  float* data = a->data;
  float* colorDataY = a->colorDataY;
  float* colorDataFY = a->colorDataFY;

  for (unsigned int i = 0; i < a->length; i++)
    {
      data[i * 3] = data[i * 3] * colorDataFY[i] / colorDataY[i];
      data[i * 3 + 1] = data[i * 3 + 1] * colorDataFY[i] / colorDataY[i];
      data[i * 3 + 2] = data[i * 3 + 2] * colorDataFY[i] / colorDataY[i];
    }

  return arg;
}

void ImageHDR::ycurve_histogram_regularization(float* colorDataY, float* colorDataFY)
{
  float yMin = colorDataY[0];
  unsigned int i = 1;
  while (yMin == 0)
    yMin = colorDataY[i++];

  for (unsigned int i = 0; i < width * height; i++)
    if (colorDataY[i] != 0 && colorDataY[i] < yMin)
      yMin = colorDataY[i];

  for (unsigned int i = 0; i < width * height; i++)
    if (colorDataY[i] == 0)
      colorDataY[i] = yMin;


  std::thread tab_t[_MT_];
  MT_histogram_regularization tab_a[_MT_];

  unsigned int id;
  unsigned int tab_length = width * height;
  unsigned int block_size = tab_length / _MT_;
  for (id = 0; id < _MT_; id++) {
    tab_a[id].data = data + (id * block_size * 3);
    tab_a[id].length = block_size;
    tab_a[id].colorDataY = colorDataY + (id * block_size);
    tab_a[id].colorDataFY = colorDataFY + (id * block_size);


    if (id == (_MT_ - 1))
      tab_a[id].length = tab_length - ((_MT_ - 1) * block_size);

    tab_t[id] = std::thread(histogram_regularization_MT, (void*)(tab_a + id));
  }

  for (id = 0; id < _MT_; id++) {
    tab_t[id].join();
  }

}

void ImageHDR::yCurve(float s, float b, float m, float w, float h)
{
  if (linear)
    {
      linear_to_non_linear();
      linear = false;
    }


  float* colorDataY = Conversion::sRGB_to_Y_of_XYZ(data, width * height);


  float yMax = colorDataY[0];
  for (unsigned int i = 0; i < width * height; i++)
    if (yMax < colorDataY[i])
      yMax = colorDataY[i];

  yMax = yMax * 100;


  YCurve yc(s, b, m, w, h, yMax);


  Eigen::MatrixXf* points = yc.evalpts(100);
  Eigen::RowVectorXf y = (*points).col(0) / 100;
  Eigen::RowVectorXf fy = (*points).col(1) / 100;
  delete(points);


  // TODO - try to optimize ?!
  // The index of the search method in utils.cpp could be calculated or determined ?
  float* colorDataFY = Utils::interp(colorDataY, width * height, y, fy);


  ycurve_histogram_regularization(colorDataY, colorDataFY);
  
  delete[](colorDataY);
  delete[](colorDataFY);
}

#else

void ImageHDR::ycurve_histogram_regularization(float* colorDataY, float* colorDataFY)
{
  float yMin = colorDataY[0];
  unsigned int i = 1;
  while (yMin == 0)
    yMin = colorDataY[i++];

  for (unsigned int i = 0; i < width * height; i++)
    if (colorDataY[i] != 0 && colorDataY[i] < yMin)
      yMin = colorDataY[i];

  for (unsigned int i = 0; i < width * height; i++)
    if (colorDataY[i] == 0)
      colorDataY[i] = yMin;

  for (unsigned int i = 0; i < width * height; i++)
    {
      data[i * 3] = data[i * 3] * colorDataFY[i] / colorDataY[i];
      data[i * 3 + 1] = data[i * 3 + 1] * colorDataFY[i] / colorDataY[i];
      data[i * 3 + 2] = data[i * 3 + 2] * colorDataFY[i] / colorDataY[i];
    }

}

void ImageHDR::yCurve(float s, float b, float m, float w, float h)
{
  if (linear)
    {
      linear_to_non_linear();
      linear = false;
    }

  float* colorDataY = Conversion::sRGB_to_Y_of_XYZ(data, width * height);


  float yMax = colorDataY[0];
  for (unsigned int i = 0; i < width * height; i++)
    if (yMax < colorDataY[i])
      yMax = colorDataY[i];

  yMax = yMax * 100;

  YCurve yc(s, b, m, w, h, yMax);


  Eigen::MatrixXf* points = yc.evalpts(100);
  Eigen::RowVectorXf y = (*points).col(0) / 100;
  Eigen::RowVectorXf fy = (*points).col(1) / 100;
  delete(points);


  float* colorDataFY = Utils::interp(colorDataY, width * height, y, fy);


  ycurve_histogram_regularization(colorDataY, colorDataFY);

  delete[](colorDataY);
  delete[](colorDataFY);
}

#endif



/****************************************/
/************** LIGHTNESSMASK ***********/
/****************************************/

#ifdef _MT_

void* lightness_MT(void* arg)
{
  MT_lightnessMask* a = (MT_lightnessMask*)arg;

  float* data = a->data;
  float* colorDataY = a->colorDataY;
  bool* mask = a->mask;

  float rangeMask[5][2] =
    {
     {0.0f, 0.2f},
     {0.2f, 0.4f},
     {0.4f, 0.6f},
     {0.6f, 0.8f},
     {0.8f, 1.0f}
    };

  unsigned int maskColor[5][3] =
    {
     {0, 0, 1},
     {0, 1, 1},
     {0, 1, 0},
     {1, 1, 0},
     {1, 0, 0}
    };

  for (unsigned int i = 0; i < a->length; i++)
    {
      for (unsigned int j = 0; j < 5; j++)
	if (mask[j])
	  if (colorDataY[i] >= rangeMask[j][0] && colorDataY[i] <= rangeMask[j][1])
	    {
	      data[i * 3] = (float)(maskColor[j][0]);
	      data[i * 3 + 1] = (float)(maskColor[j][1]);
	      data[i * 3 + 2] = (float)(maskColor[j][2]);
	    }
    }

  return arg;
}

void ImageHDR::lightnessMask(bool s, bool b, bool m, bool w, bool h)
{
  bool mask[5] = { s, b, m, w, h };

  if (linear)
    {
      linear_to_non_linear();
      linear = false;
    }

  float* colorDataY = Conversion::sRGB_to_Y_of_XYZ(data, width * height);

  std::thread tab_t[_MT_];
  MT_lightnessMask tab_a[_MT_];

  unsigned int id;
  unsigned int tab_length = width * height;
  unsigned int block_size = tab_length / _MT_;
  for (id = 0; id < _MT_; id++) {
    tab_a[id].data = data + (id * block_size * 3);
    tab_a[id].length = block_size;
    tab_a[id].colorDataY = colorDataY + (id * block_size);
    tab_a[id].mask = mask;


    if (id == (_MT_ - 1))
      tab_a[id].length = tab_length - ((_MT_ - 1) * block_size);

    tab_t[id] = std::thread(lightness_MT, (void*)(tab_a + id));
  }

  for (id = 0; id < _MT_; id++) {
    tab_t[id].join();
  }

  delete[](colorDataY);
}

#else

void ImageHDR::lightnessMask(bool s, bool b, bool m, bool w, bool h)
{
  bool mask[5] = { s, b, m, w, h };

  float rangeMask[5][2] =
    {
     {0.0f, 0.2f},
     {0.2f, 0.4f},
     {0.4f, 0.6f},
     {0.6f, 0.8f},
     {0.8f, 1.0f}
    };

  unsigned int maskColor[5][3] =
    {
     {0, 0, 1},
     {0, 1, 1},
     {0, 1, 0},
     {1, 1, 0},
     {1, 0, 0}
    };

  if (linear)
    {
      linear_to_non_linear();
      linear = false;
    }

  float* colorDataY = Conversion::sRGB_to_Y_of_XYZ(data, width * height);

  for (unsigned int i = 0; i < width * height; i++)
    {
      for (unsigned int j = 0; j < 5; j++)
	if (mask[j])
	  if (colorDataY[i] >= rangeMask[j][0] && colorDataY[i] <= rangeMask[j][1])
	    {
	      data[i * 3] = (float)(maskColor[j][0]);
	      data[i * 3 + 1] = (float)(maskColor[j][1]);
	      data[i * 3 + 2] = (float)(maskColor[j][2]);
	    }
    }

  delete[](colorDataY);
}

#endif



/****************************************/
/************** SATURATION **************/
/****************************************/

#ifdef _MT_

void* saturation_MT(void* arg)
{
  MT_saturation* a = (MT_saturation*)arg;

  float* dataLab = a->dataLab;

  for (unsigned int i = 0; i < a->length; i++)
    {
      float a_of_Lab = dataLab[i * 3 + 1];
      float b_of_Lab = dataLab[i * 3 + 2];
      dataLab[i * 3 + 1] = Conversion::Lab_to_C_of_LCH(a_of_Lab, b_of_Lab);
      // Application de la saturation
      dataLab[i * 3 + 1] = powf(dataLab[i * 3 + 1] / 100.0f, a->gamma) * 100.0f;
      dataLab[i * 3 + 2] = Conversion::Lab_to_H_of_LCH(a_of_Lab, b_of_Lab);
    }

  return arg;
}

void ImageHDR::saturation(float s)
{
  float gamma = 1.0f / ((s / 25.0f) + 1.0f);
  if (s < 0)
    gamma = (-s / 25.0f) + 1.0f;

  if (!linear)
    {
      non_linear_to_linear();
      linear = false;
    }

  float* dataLab = Conversion::sRGB_to_Lab(data, width * height);

  std::thread tab_t[_MT_];
  MT_saturation tab_a[_MT_];

  unsigned int id;
  unsigned int tab_length = width * height;
  unsigned int block_size = tab_length / _MT_;
  for (id = 0; id < _MT_; id++) {
    tab_a[id].dataLab = dataLab + (id * block_size * 3);
    tab_a[id].length = block_size;
    tab_a[id].gamma = gamma;

    if (id == (_MT_ - 1))
      tab_a[id].length = tab_length - ((_MT_ - 1) * block_size);

    tab_t[id] = std::thread(saturation_MT, (void*)(tab_a + id));
  }

  for (id = 0; id < _MT_; id++) {
    tab_t[id].join();
  }


  delete[](data);
  data = dataLab;

  linear = false;
  colorspace = Colorspace::LCH;
}

#else

void ImageHDR::saturation(float s)
{
  float gamma = 1.0f / ((s / 25.0f) + 1.0f);
  if (s < 0)
    gamma = (-s / 25.0f) + 1.0f;

  if (!linear)
    {
      non_linear_to_linear();
      linear = false;
    }

  float* dataLab = Conversion::sRGB_to_Lab(data, width * height);

  for (unsigned int i = 0; i < width * height; i++)
    {
      float a = dataLab[i * 3 + 1];
      float b = dataLab[i * 3 + 2];
      dataLab[i * 3 + 1] = Conversion::Lab_to_C_of_LCH(a, b);
      // Application de la saturation
      dataLab[i * 3 + 1] = powf(dataLab[i * 3 + 1] / 100.0f, gamma) * 100.0f;
      dataLab[i * 3 + 2] = Conversion::Lab_to_H_of_LCH(a, b);
    }

  delete[](data);
  data = dataLab;

  linear = false;
  colorspace = Colorspace::LCH;
}

#endif



/*************************************/
/************ COLOREDITOR ************/
/*************************************/

#ifdef _MT_

void* colorEditor_MT(void* arg)
{

  MT_colorEditor* a = (MT_colorEditor*)arg;

  float* data = a->data;
  unsigned int length = a->length;

  unsigned int colorspace = a->colorspace;
  bool linear = a->linear;

  float lMin = a->lMin, lMax = a->lMax;
  float cMin = a->cMin, cMax = a->cMax;
  float hMin = a->hMin, hMax = a->hMax;
  float tolerance = a->tolerance;
  float edit_hue = a->edit_hue;
  float edit_exposure = a->edit_exposure;
  float edit_contrast = a->edit_contrast;
  float edit_saturation = a->edit_saturation;
  float hueTolerance = tolerance * 360.0f;
  float chromaTolerance = tolerance * 100.0f;
  float lightTolerance = tolerance * 100.0f;
  bool mask = a->mask;


  float* dataLCH = NULL;
  float* minMask = NULL;
  float* compMask = NULL;


  // not the default parameter
  if (!(lMin == 0.0f && lMax == 100.0f
	&& cMin == 0.0f && cMax == 100.0f
	&& hMin == 0.0f && hMax == 360.0f
	&& tolerance == 0.1f
	&& edit_hue == 0.0f
	&& edit_exposure == 0.0f
	&& edit_contrast == 0.0f
	&& edit_saturation == 0.0f
	&& mask == false))
    {
      if (colorspace == Colorspace::RGB)
	{
	  if (!linear)
	    {
	      data = Conversion::non_linear_to_linear(data, length * 3);
	      linear = true;
	    }

	  dataLCH = Conversion::sRGB_to_Lab(data, length);

	  for (unsigned int i = 0; i < length; i++)
	    {
	      float a = dataLCH[i * 3 + 1];
	      float b = dataLCH[i * 3 + 2];
	      dataLCH[i * 3 + 1] = Conversion::Lab_to_C_of_LCH(a, b);
	      dataLCH[i * 3 + 2] = Conversion::Lab_to_H_of_LCH(a, b);
	    }
	}
      else
	dataLCH = data;

      float* lChannel = new float[length];
      float* cChannel = new float[length];
      float* hChannel = new float[length];

      for (unsigned int i = 0; i < length; i++)
	{
	  lChannel[i] = dataLCH[i * 3];
	  cChannel[i] = dataLCH[i * 3 + 1];
	  hChannel[i] = dataLCH[i * 3 + 2];
	}

      // Récupération du max du canal L et C
      float lMaxChannel = dataLCH[0];
      float cMaxChannel = dataLCH[1];
      for (unsigned int i = 1; i < length; i++)
	{
	  if (dataLCH[i * 3] > lMaxChannel)
	    lMaxChannel = dataLCH[i * 3];
	  if (dataLCH[i * 3 + 1] > cMaxChannel)
	    cMaxChannel = dataLCH[i * 3 + 1];
	}

      if (lMaxChannel < 100.0f)
	lMaxChannel = 100.0f;
      if (cMaxChannel < 100.0f)
	cMaxChannel = 100.0f;

      lMax = lMax * lMaxChannel / 100.0f;
      cMax = cMax * cMaxChannel / 100.0f;

      float* lightnessMask = Utils::NPlinearWeightMask(lChannel, length, lMin, lMax, lightTolerance);

      float* chromaMask = Utils::NPlinearWeightMask(cChannel, length, cMin, cMax, chromaTolerance);


      float* hueMask = Utils::NPlinearWeightMask(hChannel, length, hMin, hMax, hueTolerance);


      minMask = new float[length];
      compMask = new float[length];


      for (unsigned int i = 0; i < length; i++)
	{
	  minMask[i] = lightnessMask[i];
	  if (chromaMask[i] < minMask[i])
	    minMask[i] = chromaMask[i];
	  if (hueMask[i] < minMask[i])
	    minMask[i] = hueMask[i];
	  compMask[i] = 1.0f - minMask[i];
	}

      delete[](lightnessMask);
      delete[](chromaMask);
      delete[](hueMask);

      float hueShift = edit_hue;

      for (unsigned int i = 0; i < length; i++) {
	float oldValue = hChannel[i];
	hChannel[i] = oldValue + hueShift;
	while (hChannel[i] < 0.0f)
	  hChannel[i] += 360.0f;
	while (hChannel[i] >= 360.0f)
	  hChannel[i] -= 360.0f;
	hChannel[i] = hChannel[i] * minMask[i] + oldValue * compMask[i];
      }


      float saturation = edit_saturation;
      float gamma = 1.0f / ((saturation / 25.0f) + 1.0f);
      if (saturation < 0)
	gamma = (-saturation / 25.0f) + 1.0f;


      for (unsigned int i = 0; i < length; i++) {
	cChannel[i] = powf(cChannel[i] / 100.0f, gamma) * 100 * minMask[i] + cChannel[i] * compMask[i];
      }



      float* colorLCH = new float[length * 3];

      for (unsigned int i = 0; i < length; i++)
	{
	  colorLCH[i * 3] = lChannel[i];
	  colorLCH[i * 3 + 1] = cChannel[i];
	  colorLCH[i * 3 + 2] = hChannel[i];
	}

      delete[](lChannel);
      delete[](cChannel);
      delete[](hChannel);

      float ev = edit_exposure;
      float* colorRGB = NULL;

      if (ev != 0)
	{
	  colorRGB = Conversion::LCH_to_sRGB(colorLCH, length);

	  float* colorRGBev = new float[length * 3];
	  float coeff = powf(2, ev);

	  for (unsigned int i = 0; i < length; i++)
	    {
	      colorRGBev[i * 3] = colorRGB[i * 3] * coeff * minMask[i];
	      colorRGBev[i * 3 + 1] = colorRGB[i * 3 + 1] * coeff * minMask[i];
	      colorRGBev[i * 3 + 2] = colorRGB[i * 3 + 2] * coeff * minMask[i];

	      colorRGB[i * 3] = colorRGB[i * 3] * compMask[i] + colorRGBev[i * 3];
	      colorRGB[i * 3 + 1] = colorRGB[i * 3 + 1] * compMask[i] + colorRGBev[i * 3 + 1];
	      colorRGB[i * 3 + 2] = colorRGB[i * 3 + 2] * compMask[i] + colorRGBev[i * 3 + 2];

	    }

	  delete[](colorRGBev);

	}

      if (edit_contrast != 0)
	{
	  float contrast = edit_contrast / 100.0f;
	  float maxContrastFactor = 2.0f;
	  float scalingFactor = (1.0f - contrast) + maxContrastFactor * contrast;
	  if (contrast < 0.0f)
	    {
	      contrast = -contrast;
	      scalingFactor = 1.0f / scalingFactor;
	    }

	  float pivot = powf(2, ev) * (lMin + lMax) / 2.0f / 100.0f;


	  if (colorRGB == NULL)
	    colorRGB = Conversion::LCH_to_sRGB(colorLCH, length);
	  
	  float* colorRGB2 = Conversion::linear_to_non_linear(colorRGB, length * 3);
	  delete[](colorRGB);
	  colorRGB=colorRGB2;


	  float* colorRGBcon = new float[length * 3];
      
	  for (unsigned int i = 0; i < length; i++)
	    {
	      colorRGBcon[i * 3] = (colorRGB[i * 3] - pivot) * scalingFactor + pivot;
	      colorRGBcon[i * 3 + 1] = (colorRGB[i * 3 + 1] - pivot) * scalingFactor + pivot;
	      colorRGBcon[i * 3 + 2] = (colorRGB[i * 3 + 2] - pivot) * scalingFactor + pivot;
	      colorRGB[i * 3] = colorRGBcon[i * 3] * minMask[i] + colorRGB[i * 3] * compMask[i];
	      colorRGB[i * 3 + 1] = colorRGBcon[i * 3 + 1] * minMask[i] + colorRGB[i * 3 + 1] * compMask[i];
	      colorRGB[i * 3 + 2] = colorRGBcon[i * 3 + 2] * minMask[i] + colorRGB[i * 3 + 2] * compMask[i];
	    }

	  delete[](colorRGBcon);

	  colorRGB2 = Conversion::non_linear_to_linear(colorRGB, length * 3);
	  delete[](colorRGB);
	  colorRGB=colorRGB2;

	}

      if (colorRGB == NULL)
	colorRGB = Conversion::LCH_to_sRGB(colorLCH, length);

      for (unsigned int i = 0; i < length * 3; i++)
	{
	  data[i] = colorRGB[i];
	}

      delete[](colorRGB);
      delete[](colorLCH);
      colorspace = Colorspace::RGB;
      linear = true;

    }
  else
    {
      //TODO - To test for memory leak ?
      if (colorspace == Colorspace::LCH)
	{
	  float* colorRGB = Conversion::LCH_to_sRGB(data, length);

	  for (unsigned int i = 0; i < length * 3; i++)
	    {
	      data[i] = colorRGB[i];
	    }

	  delete[](colorRGB);

	  colorspace = Colorspace::RGB;
	  linear = true;
	}
    }

  if (mask)
    {

      for (unsigned int i = 0; i < length; i++)
	{
	  data[i * 3] = minMask[i];
	  data[i * 3 + 1] = minMask[i];
	  data[i * 3 + 2] = minMask[i];
	}
      colorspace = Colorspace::RGB;
      linear = false;
    }

  delete[](minMask);
  delete[](compMask);

  return arg;

}

void ImageHDR::colorEditor(float* selection_lightness, float* selection_chroma, float* selection_hue, float tolerance, float edit_hue, float edit_exposure, float edit_contrast, float edit_saturation, bool mask)
{
  float lMin = selection_lightness[0], lMax = selection_lightness[1];
  float cMin = selection_chroma[0], cMax = selection_chroma[1];
  float hMin = selection_hue[0], hMax = selection_hue[1];

  std::thread tab_t[_MT_];
  MT_colorEditor tab_a[_MT_];

  unsigned int id;
  unsigned int length = width * height;
  unsigned int block_size = length / _MT_;
  for (id = 0; id < _MT_; id++) {
    tab_a[id].data = data + (id * block_size * 3);
    tab_a[id].length = block_size;
    tab_a[id].colorspace = colorspace;
    tab_a[id].linear = linear;

    tab_a[id].lMin = lMin;
    tab_a[id].lMax = lMax;
    tab_a[id].cMin = cMin;
    tab_a[id].cMax = cMax;
    tab_a[id].hMin = hMin;
    tab_a[id].hMax = hMax;
    tab_a[id].tolerance = tolerance;
    tab_a[id].edit_hue = edit_hue;
    tab_a[id].edit_exposure = edit_exposure;
    tab_a[id].edit_contrast = edit_contrast;
    tab_a[id].edit_saturation = edit_saturation;
    tab_a[id].mask = mask;

    if (id == (_MT_ - 1))
      tab_a[id].length = length - ((_MT_ - 1) * block_size);

    tab_t[id] = std::thread(colorEditor_MT, (void*)(tab_a + id));
  }

  for (id = 0; id < _MT_; id++) {
    tab_t[id].join();
  }

  colorspace = Colorspace::RGB;
  linear = true;

}

#else

void ImageHDR::colorEditor(float* selection_lightness, float* selection_chroma, float* selection_hue, float tolerance, float edit_hue, float edit_exposure, float edit_contrast, float edit_saturation, bool mask)
{
  float lMin = selection_lightness[0], lMax = selection_lightness[1];
  float cMin = selection_chroma[0], cMax = selection_chroma[1];
  float hMin = selection_hue[0], hMax = selection_hue[1];
  float hueTolerance = tolerance * 360.0f;
  float chromaTolerance = tolerance * 100.0f;
  float lightTolerance = tolerance * 100.0f;
  float* dataLCH = NULL;
  float* minMask = NULL;
  float* compMask = NULL;


  // not the default parameter
  if (!(selection_lightness[0] == 0.0f && selection_lightness[1] == 100.0f
	&& selection_chroma[0] == 0.0f && selection_chroma[1] == 100.0f
	&& selection_hue[0] == 0.0f && selection_hue[1] == 360.0f
	&& tolerance == 0.1f
	&& edit_hue == 0.0f
	&& edit_exposure == 0.0f
	&& edit_contrast == 0.0f
	&& edit_saturation == 0.0f
	&& mask == false))
    {
      if (colorspace == Colorspace::RGB)
	{
	  if (!linear)
	    {
	      non_linear_to_linear();
	      linear = true;
	    }

	  dataLCH = Conversion::sRGB_to_Lab(data, width * height);

	  for (unsigned int i = 0; i < width * height; i++)
	    {
	      float a = dataLCH[i * 3 + 1];
	      float b = dataLCH[i * 3 + 2];
	      dataLCH[i * 3 + 1] = Conversion::Lab_to_C_of_LCH(a, b);
	      dataLCH[i * 3 + 2] = Conversion::Lab_to_H_of_LCH(a, b);
	    }
	}
      else
	dataLCH = data;

      float* lChannel = new float[width * height];
      float* cChannel = new float[width * height];
      float* hChannel = new float[width * height];
      for (unsigned int i = 0; i < width * height; i++)
	{
	  lChannel[i] = dataLCH[i * 3];
	  cChannel[i] = dataLCH[i * 3 + 1];
	  hChannel[i] = dataLCH[i * 3 + 2];
	}

      // Récupération du max du canal L et C
      float lMaxChannel = dataLCH[0];
      float cMaxChannel = dataLCH[1];
      for (unsigned int i = 1; i < width * height; i++)
	{
	  if (dataLCH[i * 3] > lMaxChannel)
	    lMaxChannel = dataLCH[i * 3];
	  if (dataLCH[i * 3 + 1] > cMaxChannel)
	    cMaxChannel = dataLCH[i * 3 + 1];
	}

      if (lMaxChannel < 100.0f)
	lMaxChannel = 100.0f;
      if (cMaxChannel < 100.0f)
	cMaxChannel = 100.0f;

      lMax = lMax * lMaxChannel / 100.0f;
      cMax = cMax * cMaxChannel / 100.0f;

      float* lightnessMask = Utils::NPlinearWeightMask(lChannel, width * height, lMin, lMax, lightTolerance);


      float* chromaMask = Utils::NPlinearWeightMask(cChannel, width * height, cMin, cMax, chromaTolerance);


      float* hueMask = Utils::NPlinearWeightMask(hChannel, width * height, hMin, hMax, hueTolerance);


      minMask = new float[width * height];
      compMask = new float[width * height];
      for (unsigned int i = 0; i < width * height; i++)
	{
	  minMask[i] = lightnessMask[i];
	  if (chromaMask[i] < minMask[i])
	    minMask[i] = chromaMask[i];
	  if (hueMask[i] < minMask[i])
	    minMask[i] = hueMask[i];
	  compMask[i] = 1.0f - minMask[i];
	}

      delete[](lightnessMask);
      delete[](chromaMask);
      delete[](hueMask);

      float hueShift = edit_hue;

      for (unsigned int i = 0; i < width * height; i++) {
	float oldValue = hChannel[i];
	hChannel[i] = oldValue + hueShift;
	while (hChannel[i] < 0.0f)
	  hChannel[i] += 360.0f;
	while (hChannel[i] >= 360.0f)
	  hChannel[i] -= 360.0f;
	hChannel[i] = hChannel[i] * minMask[i] + oldValue * compMask[i];
      }


      float saturation = edit_saturation;
      float gamma = 1.0f / ((saturation / 25.0f) + 1.0f);
      if (saturation < 0)
	gamma = (-saturation / 25.0f) + 1.0f;

      for (unsigned int i = 0; i < width * height; i++) {
	cChannel[i] = powf(cChannel[i] / 100.0f, gamma) * 100 * minMask[i] + cChannel[i] * compMask[i];
      }



      float* colorLCH = new float[width * height * 3];
      for (unsigned int i = 0; i < width * height; i++)
	{
	  colorLCH[i * 3] = lChannel[i];
	  colorLCH[i * 3 + 1] = cChannel[i];
	  colorLCH[i * 3 + 2] = hChannel[i];
	}

      delete[](lChannel);
      delete[](cChannel);
      delete[](hChannel);

      float ev = edit_exposure;
      float* colorRGB = NULL;

      if (ev != 0)
	{
	  colorRGB = Conversion::LCH_to_sRGB(colorLCH, width * height);

	  float* colorRGBev = new float[width * height * 3];
	  float coeff = powf(2, ev);
	  for (unsigned int i = 0; i < width * height; i++)
	    {
	      colorRGBev[i * 3] = colorRGB[i * 3] * coeff * minMask[i];
	      colorRGBev[i * 3 + 1] = colorRGB[i * 3 + 1] * coeff * minMask[i];
	      colorRGBev[i * 3 + 2] = colorRGB[i * 3 + 2] * coeff * minMask[i];

	      colorRGB[i * 3] = colorRGB[i * 3] * compMask[i] + colorRGBev[i * 3];
	      colorRGB[i * 3 + 1] = colorRGB[i * 3 + 1] * compMask[i] + colorRGBev[i * 3 + 1];
	      colorRGB[i * 3 + 2] = colorRGB[i * 3 + 2] * compMask[i] + colorRGBev[i * 3 + 2];

	    }

	  delete[](colorRGBev);

	}

      if (edit_contrast != 0)
	{
	  float contrast = edit_contrast / 100.0f;
	  float maxContrastFactor = 2.0f;
	  float scalingFactor = (1.0f - contrast) + maxContrastFactor * contrast;
	  if (contrast < 0.0f)
	    {
	      contrast = -contrast;
	      scalingFactor = 1.0f / scalingFactor;
	    }

	  float pivot = powf(2, ev) * (lMin + lMax) / 2.0f / 100.0f;


	  if (colorRGB == NULL)
	    colorRGB = Conversion::LCH_to_sRGB(colorLCH, width * height);
	  
	  float* colorRGB2 = Conversion::linear_to_non_linear(colorRGB, width * height * 3);
	  delete[](colorRGB);
	  colorRGB=colorRGB2;


	  float* colorRGBcon = new float[width * height * 3];
	  for (unsigned int i = 0; i < width * height; i++)
	    {
	      colorRGBcon[i * 3] = (colorRGB[i * 3] - pivot) * scalingFactor + pivot;
	      colorRGBcon[i * 3 + 1] = (colorRGB[i * 3 + 1] - pivot) * scalingFactor + pivot;
	      colorRGBcon[i * 3 + 2] = (colorRGB[i * 3 + 2] - pivot) * scalingFactor + pivot;
	      colorRGB[i * 3] = colorRGBcon[i * 3] * minMask[i] + colorRGB[i * 3] * compMask[i];
	      colorRGB[i * 3 + 1] = colorRGBcon[i * 3 + 1] * minMask[i] + colorRGB[i * 3 + 1] * compMask[i];
	      colorRGB[i * 3 + 2] = colorRGBcon[i * 3 + 2] * minMask[i] + colorRGB[i * 3 + 2] * compMask[i];
	    }

	  delete[](colorRGBcon);

	  colorRGB2 = Conversion::non_linear_to_linear(colorRGB, width * height * 3);
	  delete[](colorRGB);
	  colorRGB=colorRGB2;

	}

      if (colorRGB == NULL)
	colorRGB = Conversion::LCH_to_sRGB(colorLCH, width * height);
      
      for (unsigned int i = 0; i < width*height * 3; i++)
	{
	  data[i] = colorRGB[i];
	}

      delete[](colorRGB);
      delete[](colorLCH);
      colorspace = Colorspace::RGB;
      linear = true;

    }
  else
    {
      if (colorspace == Colorspace::LCH)
	{
	  float* colorRGB = Conversion::LCH_to_sRGB(data, width * height);
	  
	  for (unsigned int i = 0; i < width*height * 3; i++)
	    {
	      data[i] = colorRGB[i];
	    }

	  delete[](colorRGB);
	  colorspace = Colorspace::RGB;
	  linear = true;
	}
    }

  if (mask)
    {
      for (unsigned int i = 0; i < width * height; i++)
	{
	  data[i * 3] = minMask[i];
	  data[i * 3 + 1] = minMask[i];
	  data[i * 3 + 2] = minMask[i];
	}
      colorspace = Colorspace::RGB;
      linear = false;
    }

  delete[](minMask);
  delete[](compMask);

}

#endif


/* Private methods */


Eigen::VectorXf ImageHDR::to_EigenVector() const
{
  Eigen::VectorXf v(width * height * 3);

  for (unsigned int i = 0; i < width; i++)
    v(i) = data[i];

  return v;
}