---

color.h

Go to the documentation of this file.

#ifndef __COLOR__
#define __COLOR__

#include <math.h>

#include <iterator>
#include <list>
#include <vector>

#ifndef WITHOUT_LIMITS
#include <limits>
#endif

#include "io_tools.h"
#include "msg_stream.h"

#ifndef WITHOUT_STATIC_CONST  // in case of trouble with gcc 2.96
#define CONST_STATIC const
#else
#define CONST_STATIC
#endif

/*

  ###############
  # typedef ... #
  ###############
  
 */

template<typename RGB_converter> class Basic_color;

template<typename Real> class RGB_color;
template<typename Real> class HSV_color;
template<typename Real> class XYZ_color;
template<typename Real> class Luv_color;
template<typename Real> class Grey_level_color;

template<typename Real> class RGB_converter_to_RGB;
template<typename Real> class RGB_converter_to_HSV;
template<typename Real> class RGB_converter_to_XYZ;
template<typename Real> class RGB_converter_to_Luv;
template<typename Real> class RGB_converter_to_grey_level;

typedef RGB_color<float>         RGB_f_color;
typedef RGB_color<double>        RGB_d_color;

typedef HSV_color<float>         HSV_f_color;
typedef HSV_color<double>        HSV_d_color;

typedef XYZ_color<float>         XYZ_f_color;
typedef XYZ_color<double>        XYZ_d_color;

typedef Luv_color<float>         Luv_f_color;
typedef Luv_color<double>        Luv_d_color;

typedef Grey_level_color<float>  Grey_level_f_color;
typedef Grey_level_color<double> Grey_level_d_color;


template<typename RGB_converter>
std::ostream& operator<<(std::ostream& s,
                         const Basic_color<RGB_converter>& c);

/*

  #####################
  # class Basic_color #
  #####################

*/


template<typename RGB_converter>
class Basic_color{

public:

  typedef typename RGB_converter::component_type component_type;
  typedef component_type                         distance_type;

  typedef RGB_converter  converter_type;
  
  Basic_color();
  
  Basic_color(const char r,  const char g,  const char b);
  Basic_color(const float r, const float g, const float b);
  Basic_color(const double r,const double g,const double b);
  
  Basic_color(const Basic_color<RGB_converter>& c);
  
  inline Basic_color<RGB_converter>&
  operator=(const Basic_color<RGB_converter>& c);

  inline bool operator==(const Basic_color<RGB_converter>& c) const;
  inline bool operator!=(const Basic_color<RGB_converter>& c) const;
  

  
  inline void get_RGB(char* const r,  char* const g,  char* const b) const;
  inline void get_RGB(float* const r, float* const g, float* const b) const;
  inline void get_RGB(double* const r,double* const g,double* const b) const;  


  
  inline void get_space_position(component_type* const x_space,
                                 component_type* const y_space,
                                 component_type* const z_space) const;

  template<typename Space_vector>
  inline void get_space_position(Space_vector* const space_vector) const;

  

  
  inline void set_space_position(const component_type x_space,
                                 const component_type y_space,
                                 const component_type z_space);

  template<typename Space_vector>
  inline void set_space_position(const Space_vector& space_vector);

  inline bool is_not_a_color() const;

  inline void write_to_bytes(std::ostream& out) const;

  inline void read_from_bytes(std::istream& in);

  inline static distance_type
  square_distance(const Basic_color<RGB_converter>& c1,
                  const Basic_color<RGB_converter>& c2);

  inline static distance_type
  distance(const Basic_color<RGB_converter>& c1,
           const Basic_color<RGB_converter>& c2);
  

  
  template<typename Basic_color_iterator>
  static void mean(Basic_color_iterator first,
                   Basic_color_iterator last,
                   Basic_color<RGB_converter>* const result);

  template<typename Basic_color_iterator,typename Weight_iterator>
  static void mean(Basic_color_iterator first_clr,
                   Basic_color_iterator last_clr,
                   Weight_iterator first_w,
                   Weight_iterator last_w,
                   Basic_color<RGB_converter>* const result);



  template<typename Basic_color_iterator>
  static void robust_mean(Basic_color_iterator first,
                          Basic_color_iterator last,
                          Basic_color<RGB_converter>* const result,
                          const unsigned int n_iter = 1);

  template<typename Basic_color_iterator,typename Weight_iterator>
  static void robust_mean(Basic_color_iterator first_clr,
                          Basic_color_iterator last_clr,
                          Weight_iterator first_w,
                          Weight_iterator last_w,
                          Basic_color<RGB_converter>* const result,
                          const unsigned int n_iter = 1);

  static void main_conversion_output(std::ostream& out = std::cout) ;


  
  template<typename Basic_color_iterator>
  inline static distance_type
  square_deviation(const Basic_color_iterator first,
                   const Basic_color_iterator last);

  template<typename Basic_color_iterator>
  static distance_type
  square_deviation(const Basic_color_iterator first,
                   const Basic_color_iterator last,
                   Basic_color<RGB_converter>* const mean_color);



  static const Basic_color<RGB_converter> black;
  static const Basic_color<RGB_converter> grey;
  static const Basic_color<RGB_converter> white;
  
  static const Basic_color<RGB_converter> red;
  static const Basic_color<RGB_converter> green;
  static const Basic_color<RGB_converter> blue;

  static const Basic_color<RGB_converter> cyan;
  static const Basic_color<RGB_converter> magenta;
  static const Basic_color<RGB_converter> yellow;

  static void cycle(const double alpha,
                    Basic_color<RGB_converter>* const result,
                    const Basic_color<RGB_converter>& clr_0 = red,
                    const Basic_color<RGB_converter>& clr_1 = blue,
                    const Basic_color<RGB_converter>& clr_2 = yellow,
                    const Basic_color<RGB_converter>& clr_3 = magenta);
  
  friend std::ostream&
  operator<<<RGB_converter>(std::ostream& s,
                            const Basic_color<RGB_converter>& c);

private:

  
  component_type x,y,z;
};


/*

  ###################
  # class RGB_color #
  ###################

 */


template<typename Real>
class RGB_color
  :public Basic_color<RGB_converter_to_RGB<Real> >{

public:
  
  RGB_color()
    :Basic_color<RGB_converter_to_RGB<Real> >(){}
  
  RGB_color(const char r,  const char g,  const char b)
    :Basic_color<RGB_converter_to_RGB<Real> >(r,g,b){}

  RGB_color(const float r, const float g, const float b)
    :Basic_color<RGB_converter_to_RGB<Real> >(r,g,b){}

  RGB_color(const double r,const double g,const double b)
    :Basic_color<RGB_converter_to_RGB<Real> >(r,g,b){}
  
  RGB_color(const Basic_color<RGB_converter_to_RGB<Real> >& c)
    :Basic_color<RGB_converter_to_RGB<Real> >(c){}

};






/*

  ###################
  # class HSV_color #
  ###################

 */


template<typename Real>
class HSV_color
  :public Basic_color<RGB_converter_to_HSV<Real> >{

public:
  
  HSV_color()
    :Basic_color<RGB_converter_to_HSV<Real> >(){}
  
  HSV_color(const char r,  const char g,  const char b)
    :Basic_color<RGB_converter_to_HSV<Real> >(r,g,b){}
  
  HSV_color(const float r, const float g, const float b)
    :Basic_color<RGB_converter_to_HSV<Real> >(r,g,b){}
  
  HSV_color(const double r,const double g,const double b)
    :Basic_color<RGB_converter_to_HSV<Real> >(r,g,b){}
  
  HSV_color(const Basic_color<RGB_converter_to_HSV<Real> >& c)
    :Basic_color<RGB_converter_to_HSV<Real> >(c){}

};



/*

  ###################
  # class XYZ_color #
  ###################

 */


template<typename Real>
class XYZ_color
  :public Basic_color<RGB_converter_to_XYZ<Real> >{

public:
  
  XYZ_color()
    :Basic_color<RGB_converter_to_XYZ<Real> >(){}
  
  XYZ_color(const char r,  const char g,  const char b)
    :Basic_color<RGB_converter_to_XYZ<Real> >(r,g,b){}
  
  XYZ_color(const float r, const float g, const float b)
    :Basic_color<RGB_converter_to_XYZ<Real> >(r,g,b){}
  
  XYZ_color(const double r,const double g,const double b)
    :Basic_color<RGB_converter_to_XYZ<Real> >(r,g,b){}
  
  XYZ_color(const Basic_color<RGB_converter_to_XYZ<Real> >& c)
    :Basic_color<RGB_converter_to_XYZ<Real> >(c){}

};





/*

  ###################
  # class Luv_color #
  ###################

 */


template<typename Real>
class Luv_color
  :public Basic_color<RGB_converter_to_Luv<Real> >{

public:
  
  Luv_color()
    :Basic_color<RGB_converter_to_Luv<Real> >(){}
  
  Luv_color(const char r,  const char g,  const char b)
    :Basic_color<RGB_converter_to_Luv<Real> >(r,g,b){}
  
  Luv_color(const float r, const float g, const float b)
    :Basic_color<RGB_converter_to_Luv<Real> >(r,g,b){}
  
  Luv_color(const double r,const double g,const double b)
    :Basic_color<RGB_converter_to_Luv<Real> >(r,g,b){}
  
  Luv_color(const Basic_color<RGB_converter_to_Luv<Real> >& c)
    :Basic_color<RGB_converter_to_Luv<Real> >(c){}

};





/*

  ##########################
  # class Grey_level_color #
  ##########################

 */


template<typename Real>
class Grey_level_color
  :public Basic_color<RGB_converter_to_grey_level<Real> >{

public:
  
  Grey_level_color()
    :Basic_color<RGB_converter_to_grey_level<Real> >(){}

  Grey_level_color(const char c)
    :Basic_color<RGB_converter_to_grey_level<Real> >(c,c,c){}

  Grey_level_color(const float f)
    :Basic_color<RGB_converter_to_grey_level<Real> >(f,f,f){}

  Grey_level_color(const double d)
    :Basic_color<RGB_converter_to_grey_level<Real> >(d,d,d){}
  
  Grey_level_color(const char r,  const char g,  const char b)
    :Basic_color<RGB_converter_to_grey_level<Real> >(r,g,b){}
  
  Grey_level_color(const float r, const float g, const float b)
    :Basic_color<RGB_converter_to_grey_level<Real> >(r,g,b){}
  
  Grey_level_color(const double r,const double g,const double b)
    :Basic_color<RGB_converter_to_grey_level<Real> >(r,g,b){}
  
  Grey_level_color(const Basic_color<RGB_converter_to_grey_level<Real> >& c)
    :Basic_color<RGB_converter_to_grey_level<Real> >(c){}

};


/*

  ##############################
  # class RGB_converter_to_HSV #
  ##############################

*/



template<typename Real>
class RGB_converter_to_HSV{
  
public:

  typedef Real component_type;
  
  inline static void from_RGB_to_cartesian(const double r,
                                           const double g,
                                           const double b,
                                           component_type* const x,
                                           component_type* const y,
                                           component_type* const z);

  inline static void from_cartesian_to_RGB(const component_type x,
                                           const component_type y,
                                           const component_type z,
                                           double* const r,
                                           double* const g,
                                           double* const b);

private:
  inline static void RGB_to_HSV(const double r,const double g,const double b,
                                double* const h,double* const s,double* const v);
  
  inline static void HSV_to_RGB(const double h,const double s,const double v,
                                double* const r,double* const g,double* const b);
};




/*

  ##############################
  # class RGB_converter_to_RGB #
  ##############################

*/



template<typename Real>
class RGB_converter_to_RGB{
  
public:

  typedef Real component_type;
  
  inline static void from_RGB_to_cartesian(const double r,
                                           const double g,
                                           const double b,
                                           component_type* const x,
                                           component_type* const y,
                                           component_type* const z);

  inline static void from_cartesian_to_RGB(const component_type x,
                                           const component_type y,
                                           const component_type z,
                                           double* const r,
                                           double* const g,
                                           double* const b);
};




template<typename Real>
class RGB_converter_to_XYZ{
  
public:

  typedef Real component_type;
  
  inline static void from_RGB_to_cartesian(const double r,
                                           const double g,
                                           const double b,
                                           component_type* const x,
                                           component_type* const y,
                                           component_type* const z);

  inline static void from_cartesian_to_RGB(const component_type x,
                                           const component_type y,
                                           const component_type z,
                                           double* const r,
                                           double* const g,
                                           double* const b);
};




template<typename Real>
class RGB_converter_to_Luv{
  
public:

  typedef Real component_type;
  
  inline static void from_RGB_to_cartesian(const double r,
                                           const double g,
                                           const double b,
                                           component_type* const x,
                                           component_type* const y,
                                           component_type* const z);

  inline static void from_cartesian_to_RGB(const component_type x,
                                           const component_type y,
                                           const component_type z,
                                           double* const r,
                                           double* const g,
                                           double* const b);
private:

  
  static const double Xn;
  static const double Yn;
  static const double Zn;


  
  static double L2Y(const double l);
  static double Y2L(const double y);
};





template<typename Real>
class RGB_converter_to_grey_level{
  
public:

  typedef Real component_type;
  
  inline static void from_RGB_to_cartesian(const double r,
                                           const double g,
                                           const double b,
                                           component_type* const x,
                                           component_type* const y,
                                           component_type* const z);

  inline static void from_cartesian_to_RGB(const component_type x,
                                           const component_type y,
                                           const component_type z,
                                           double* const r,
                                           double* const g,
                                           double* const b);

  inline static void set_min_value(const double value);
  inline static void set_max_value(const double value);
};




/*

  #######################
  # class Color_to_grey #
  #######################

 */


template<typename Color,typename Grey>
class Color_to_grey{
public:
  Color_to_grey(){}

  Grey operator()(const Color& c){
    Grey r,g,b;
    c.get_RGB(&r,&g,&b);
    return (r+g+b)/3;
  }
};





/*

  #######################
  # class Grey_to_color #
  #######################

*/

template<typename Grey,typename Color>
class Grey_to_color{
public:
  Grey_to_color(){}

  inline Color operator()(const Grey g){ return Color(g,g,g); }
};




/*

  ########################
  # class Black_or_white #
  ########################
  
*/

template<typename Color>
class Black_or_white{
public:
  Black_or_white(){}

  inline Color operator()(const bool b){ return (b?(Color(1.0,1.0,1.0)):
                                                 (Color(0.0,0.0,0.0))); }
};









/*
  
  #############################################
  #############################################
  #############################################
  ######                                 ######
  ######   I M P L E M E N T A T I O N   ######
  ######                                 ######
  #############################################
  #############################################
  #############################################
  
*/







/*

  #####################
  # class Basic_color #
  #####################

*/

#ifndef WITHOUT_LIMITS

template<typename RGB_converter>
Basic_color<RGB_converter>::
Basic_color()
  :x(std::numeric_limits<component_type>::signaling_NaN()),
   y(std::numeric_limits<component_type>::signaling_NaN()),
   z(std::numeric_limits<component_type>::signaling_NaN()){}

#else

template<typename RGB_converter>
Basic_color<RGB_converter>::
Basic_color()
  :x(),y(),z(){

}

#endif

template<typename RGB_converter>
Basic_color<RGB_converter>::
Basic_color(const char r,const char g,const char b){

  RGB_converter::from_RGB_to_cartesian(static_cast<double>(r)/255,
                                       static_cast<double>(g)/255,
                                       static_cast<double>(b)/255,
                                       &x,&y,&z);
}


template<typename RGB_converter>
Basic_color<RGB_converter>::
Basic_color(const float r,const float g,const float b){

  RGB_converter::from_RGB_to_cartesian(static_cast<double>(r),
                                       static_cast<double>(g),
                                       static_cast<double>(b),
                                       &x,&y,&z);
}


template<typename RGB_converter>
Basic_color<RGB_converter>::
Basic_color(const double r,const double g,const double b){

  RGB_converter::from_RGB_to_cartesian(r,g,b,&x,&y,&z);
}


template<typename RGB_converter>
Basic_color<RGB_converter>::
Basic_color(const Basic_color<RGB_converter>& c)
  :x(c.x),y(c.y),z(c.z){
}


template<typename RGB_converter>
Basic_color<RGB_converter>&
Basic_color<RGB_converter>::
operator=(const Basic_color<RGB_converter>& c){

  if (this!=&c){
    x = c.x;
    y = c.y;
    z = c.z;
  }
  
  return *this;
}


template<typename RGB_converter>
bool
Basic_color<RGB_converter>::
operator==(const Basic_color<RGB_converter>& c) const{
  
  return ((x==c.x)&&(y==c.y)&&(z==c.z));
}


template<typename RGB_converter>
bool
Basic_color<RGB_converter>::
operator!=(const Basic_color<RGB_converter>& c) const{
  
  return ((x!=c.x)||(y!=c.y)||(z!=c.z));
}


template<typename RGB_converter>
void
Basic_color<RGB_converter>::
get_RGB(char* const r,char* const g,char* const b) const{

  double R,G,B;
  
  RGB_converter::from_cartesian_to_RGB(x,y,z,&R,&G,&B);

  *r = static_cast<char>(rint(R*255));
  *g = static_cast<char>(rint(G*255));
  *b = static_cast<char>(rint(B*255));  
}


template<typename RGB_converter>
void
Basic_color<RGB_converter>::
get_RGB(float* const r,float* const g,float* const b) const{

  double R,G,B;
  
  RGB_converter::from_cartesian_to_RGB(x,y,z,&R,&G,&B);

  *r = static_cast<float>(R);
  *g = static_cast<float>(G);
  *b = static_cast<float>(B);  
}


template<typename RGB_converter>
void
Basic_color<RGB_converter>::
get_RGB(double* const r,double* const g,double* const b) const{
  RGB_converter::from_cartesian_to_RGB(x,y,z,r,g,b);
}


template<typename RGB_converter>
void
Basic_color<RGB_converter>::
get_space_position(component_type* const x_space,
                   component_type* const y_space,
                   component_type* const z_space) const{
  *x_space = x;
  *y_space = y;
  *z_space = z;
}


template<typename RGB_converter>
template<typename Space_vector>
void
Basic_color<RGB_converter>::
get_space_position(Space_vector* const space_vector) const{
  
  typedef typename Space_vector::value_type type;
  
  (*space_vector)[0] = static_cast<type>(x);
  (*space_vector)[1] = static_cast<type>(y);
  (*space_vector)[2] = static_cast<type>(z);
}


template<typename RGB_converter>
void
Basic_color<RGB_converter>::
set_space_position(const component_type x_space,
                   const component_type y_space,
                   const component_type z_space){
  x = x_space;
  y = y_space;
  z = z_space;
}


template<typename RGB_converter>
template<typename Space_vector>
void
Basic_color<RGB_converter>::
set_space_position(const Space_vector& space_vector){
  x = static_cast<component_type>(space_vector[0]);
  y = static_cast<component_type>(space_vector[1]);
  z = static_cast<component_type>(space_vector[2]);
}


template<typename RGB_converter>
bool
Basic_color<RGB_converter>::
is_not_a_color() const{
  return (is_signaling_NaN(x)||is_signaling_NaN(y)||is_signaling_NaN(z));  
}


template<typename RGB_converter>
void
Basic_color<RGB_converter>::
write_to_bytes(std::ostream& out) const{

  using namespace IO_tools;
  
  int size = sizeof(component_type);

  switch (size){
    
  case 2:
    write_2_bytes(x,out); 
    write_2_bytes(y,out);
    write_2_bytes(z,out);
    return;

  case 4:
    write_4_bytes(x,out); 
    write_4_bytes(y,out);
    write_4_bytes(z,out);
    return;

  default:
    Message::error<<"Basic_color::write_to_bytes : sizeof(component_type) = "<<size
                  <<" non implemented"<<Message::done;
  }
  
  return;
}


template<typename RGB_converter>
void
Basic_color<RGB_converter>::
read_from_bytes(std::istream& in){

  using namespace IO_tools;

  int size = sizeof(component_type);

  switch (size){
    
  case 2:
    read_2_bytes(&x,in); 
    read_2_bytes(&y,in);
    read_2_bytes(&z,in);
    return;

  case 4:
    read_4_bytes(&x,in); 
    read_4_bytes(&y,in);
    read_4_bytes(&z,in);
    return;

  default:
    Message::error<<"Basic_color::read_to_bytes : sizeof(component_type) = "<<size
                  <<" non implemenrted"<<Message::done;
  }
}


template<typename RGB_converter>
typename Basic_color<RGB_converter>::distance_type
Basic_color<RGB_converter>::
square_distance(const Basic_color<RGB_converter>& c1,
                const Basic_color<RGB_converter>& c2){

  const distance_type dx = c1.x - c2.x;
  const distance_type dy = c1.y - c2.y;
  const distance_type dz = c1.z - c2.z;

  return (dx*dx + dy*dy + dz*dz);
}


template<typename RGB_converter>
typename Basic_color<RGB_converter>::distance_type
Basic_color<RGB_converter>::
distance(const Basic_color<RGB_converter>& c1,
         const Basic_color<RGB_converter>& c2){
  return sqrtf(square_distance(c1,c2));
}


template<typename RGB_converter>
template<typename Basic_color_iterator>
void
Basic_color<RGB_converter>::
mean(Basic_color_iterator first,
     Basic_color_iterator last,
     Basic_color<RGB_converter>* const result){

  unsigned int n     = 0;
  result->x = 0;
  result->y = 0;
  result->z = 0;  
  
  for(Basic_color_iterator c=first;c!=last;c++){
    result->x += c->x;
    result->y += c->y;
    result->z += c->z;
    n++;
  }

  result->x /= n;
  result->y /= n;
  result->z /= n;
}


template<typename RGB_converter>
template<typename Basic_color_iterator,typename Weight_iterator>
void
Basic_color<RGB_converter>::
mean(Basic_color_iterator first_clr,
     Basic_color_iterator last_clr,
     Weight_iterator      first_w,
     Weight_iterator      last_w,
     Basic_color<RGB_converter>* const result){

  component_type n = 0;
  result->x        = 0;
  result->y        = 0;
  result->z        = 0;
  
  Basic_color_iterator c;
  Weight_iterator w;
  for(c=first_clr,w=first_w;
      (c!=last_clr)&&(w!=last_w);
      c++,w++){
    result->x += c->x * (*w);
    result->y += c->y * (*w);
    result->z += c->z * (*w);
    n         += (*w);
  }

  if (n==0){
    Message::warning<<"Basic_color::mean : the weight sum is null\n"<<Message::done;
    *result = Basic_color<RGB_converter>();
  }
  else{
    result->x /= n;
    result->y /= n;
    result->z /= n;
  }
}


template<typename RGB_converter>
template<typename Basic_color_iterator>
void
Basic_color<RGB_converter>::
robust_mean(Basic_color_iterator first,
            Basic_color_iterator last,
            Basic_color<RGB_converter>* const result,
            const unsigned int n_iter){
  
  mean(first,last,result);

  if (n_iter==0){
    return;
  }
  
  int n             = 0;
  distance_type sum = 0;

  std::list<distance_type> dist;

  // Compute the variamce.
  
  for(Basic_color_iterator c=first;c!=last;c++){
    distance_type sq = square_distance(*c,*result);
    dist.push_back(sq);
    sum += sq;
    n++;
  }

  distance_type sq_dev;

  if (n>0) {
    sq_dev = (sum/n);
  }

#ifndef WITHOUT_LIMITS
  else{
    sq_dev = std::numeric_limits<distance_type>::signaling_NaN();
    Message::warning<<"no color to compute the variance"<<Message::done;
  }
#endif

  // Keep only the nearest ones.
  typedef typename std::iterator_traits<Basic_color_iterator>::value_type
    color_type;
  
  std::list<color_type> near_one;

  n = 0;
  Basic_color_iterator               c;
  typename std::list<distance_type>::iterator d;
  for(c=first,d=dist.begin();c!=last;c++,d++){
    if ((*d)<=sq_dev){
      near_one.push_back(*c);
      n++;
    }
  }

  
  robust_mean(near_one.begin(),near_one.end(),result,n_iter-1);
}


template<typename RGB_converter>
template<typename Basic_color_iterator,typename Weight_iterator>
void
Basic_color<RGB_converter>::
robust_mean(Basic_color_iterator first_clr,
                 Basic_color_iterator last_clr,
                 Weight_iterator first_w,
                 Weight_iterator last_w,
                 Basic_color<RGB_converter>* const result,
                 const unsigned int n_iter){

  mean(first_clr,last_clr,first_w,last_w,result);

  if (n_iter==0){
    return;
  }
  
  component_type n   = 0;
  distance_type  sum = 0;

  std::list<distance_type> dist;

  // Compute the variance.
  Basic_color_iterator c;
  Weight_iterator w;
  for(c=first_clr,w=first_w;
      (c!=last_clr)&&(w!=last_w);
      c++,w++){

    distance_type sq = square_distance(*c,*result);
    dist.push_back(sq);
    sum += (*w) * sq;
    n += *w;
  }

  distance_type sq_dev;

  if (n>0) {
    sq_dev = (sum/n);
  }

#ifndef WITHOUT_LIMITS
  else{
    sq_dev = std::numeric_limits<distance_type>::signaling_NaN();
    Message::warning<<"no color to compute the variance"<<Message::done;
  }
#endif

  // Keep only the nearest one.
  typedef typename std::iterator_traits<Basic_color_iterator>::value_type
    color_type;
  typedef typename std::iterator_traits<Weight_iterator>::value_type
    weight_type;
  
  std::list<color_type> near_clr;
  std::list<weight_type> near_w;

  n = 0;
  typename std::list<distance_type>::iterator d;
  for(c=first_clr, w=first_w, d=dist.begin();
      (c != last_clr) && (w != last_w);
      c++, w++, d++){
    
    if ((*d)<=sq_dev){
      near_clr.push_back(*c);
      near_w.push_back(*w);
      n++;
    }
  }

  
  robust_mean(near_clr.begin(),near_clr.end(),
              near_w.begin(),near_w.end(),
              result,n_iter-1);

}


template<typename RGB_converter>
template<typename Basic_color_iterator>
typename Basic_color<RGB_converter>::distance_type
Basic_color<RGB_converter>::
square_deviation(Basic_color_iterator first,
                 Basic_color_iterator last){

  Basic_color<RGB_converter> mean_color;

  return square_deviation(first,last,&mean_color);
}


template<typename RGB_converter>
template<typename Basic_color_iterator>
typename Basic_color<RGB_converter>::distance_type
Basic_color<RGB_converter>::
square_deviation(Basic_color_iterator first,
                 Basic_color_iterator last,
                 Basic_color<RGB_converter>* const mean_color){

  mean(first,last,mean_color);

  int n             = 0;
  distance_type sum = 0;
  
  for(Basic_color_iterator c=first;c!=last;c++){
    sum += square_distance(*c,*mean_color);
    n++;
  }

#ifndef WITHOUT_LIMITS
  return ((n>0)?(sum/n):(std::numeric_limits<distance_type>::signaling_NaN()));
#else
  return (sum/n);
#endif
}


template<typename RGB_converter>
void
Basic_color<RGB_converter>::
main_conversion_output(std::ostream& out){

  Basic_color<RGB_converter> clr;
  component_type r,g,b;

  clr = white;
  clr.get_RGB(&r,&g,&b);
  out<<"r: "<<r<<"\tg: "<<g<<"\tb: "<<b<<"\n"
     <<"->\tx: "<<clr.x<<"\ty: "<<clr.y<<"\tz: "<<clr.z<<"\n"<<std::endl;

  clr = grey;
  clr.get_RGB(&r,&g,&b);
  out<<"r: "<<r<<"\tg: "<<g<<"\tb: "<<b<<"\n"
     <<"->\tx: "<<clr.x<<"\ty: "<<clr.y<<"\tz: "<<clr.z<<"\n"<<std::endl;

  clr = black;
  clr.get_RGB(&r,&g,&b);
  out<<"r: "<<r<<"\tg: "<<g<<"\tb: "<<b<<"\n"
     <<"->\tx: "<<clr.x<<"\ty: "<<clr.y<<"\tz: "<<clr.z<<"\n"<<std::endl;

  clr = red;
  clr.get_RGB(&r,&g,&b);
  out<<"r: "<<r<<"\tg: "<<g<<"\tb: "<<b<<"\n"
     <<"->\tx: "<<clr.x<<"\ty: "<<clr.y<<"\tz: "<<clr.z<<"\n"<<std::endl;

  clr = green;
  clr.get_RGB(&r,&g,&b);
  out<<"r: "<<r<<"\tg: "<<g<<"\tb: "<<b<<"\n"
     <<"->\tx: "<<clr.x<<"\ty: "<<clr.y<<"\tz: "<<clr.z<<"\n"<<std::endl;
  
  clr = blue;
  clr.get_RGB(&r,&g,&b);
  out<<"r: "<<r<<"\tg: "<<g<<"\tb: "<<b<<"\n"
     <<"->\tx: "<<clr.x<<"\ty: "<<clr.y<<"\tz: "<<clr.z<<"\n"<<std::endl;
  
  clr = cyan;
  clr.get_RGB(&r,&g,&b);
  out<<"r: "<<r<<"\tg: "<<g<<"\tb: "<<b<<"\n"
     <<"->\tx: "<<clr.x<<"\ty: "<<clr.y<<"\tz: "<<clr.z<<"\n"<<std::endl;

  clr = magenta;
  clr.get_RGB(&r,&g,&b);
  out<<"r: "<<r<<"\tg: "<<g<<"\tb: "<<b<<"\n"
     <<"->\tx: "<<clr.x<<"\ty: "<<clr.y<<"\tz: "<<clr.z<<"\n"<<std::endl;
  
  clr = yellow;
  clr.get_RGB(&r,&g,&b);
  out<<"r: "<<r<<"\tg: "<<g<<"\tb: "<<b<<"\n"
     <<"->\tx: "<<clr.x<<"\ty: "<<clr.y<<"\tz: "<<clr.z<<"\n"<<std::endl;
}


template<typename RGB_converter>
const Basic_color<RGB_converter>
Basic_color<RGB_converter>::white(1.0f, 1.0f, 1.0f);

template<typename RGB_converter>
const Basic_color<RGB_converter>
Basic_color<RGB_converter>::grey(0.5f, 0.5f, 0.5f);

template<typename RGB_converter>
const Basic_color<RGB_converter>
Basic_color<RGB_converter>::black(0.0f, 0.0f, 0.0f);

template<typename RGB_converter>
const Basic_color<RGB_converter>
Basic_color<RGB_converter>::red(1.0f, 0.0f, 0.0f);

template<typename RGB_converter>
const Basic_color<RGB_converter>
Basic_color<RGB_converter>::green(0.0f, 1.0f, 0.0f);

template<typename RGB_converter>
const Basic_color<RGB_converter>
Basic_color<RGB_converter>::blue(0.0f, 0.0f, 1.0f);

template<typename RGB_converter>
const Basic_color<RGB_converter>
Basic_color<RGB_converter>::cyan(0.0f, 1.0f, 1.0f);

template<typename RGB_converter>
const Basic_color<RGB_converter>
Basic_color<RGB_converter>::magenta(1.0f, 0.0f, 1.0f);

template<typename RGB_converter>
const Basic_color<RGB_converter>
Basic_color<RGB_converter>::yellow(1.0f, 1.0f, 0.0f);


template<typename RGB_converter>
void
Basic_color<RGB_converter>::
cycle(const double alpha,
      Basic_color<RGB_converter>* const result,
      const Basic_color<RGB_converter>& clr_0,
      const Basic_color<RGB_converter>& clr_1,
      const Basic_color<RGB_converter>& clr_2,
      const Basic_color<RGB_converter>& clr_3){

  const double a = alpha - floor(alpha/(2.0*M_PI))*2.0*M_PI;
  
  std::vector<Basic_color<RGB_converter> > color(2);
  std::vector<double>                      weight(2);


  if ((a<0) || (a>2.0*M_PI)){
    Message::error<<"Basic_color<RGB_converter>::cycle : valeur de a"
                  <<" non valable : "<<a<<Message::done;
  }
  
  if ((a>=0) && (a<=M_PI_2)){
    color[0]  = clr_0;
    weight[0] = M_PI_2 - a;
    
    color[1]  = clr_1;
    weight[1] = M_PI_2 - weight[0]; 
  }

  if ((a>M_PI_2) && (a<=M_PI)){
    color[0]  = clr_1;
    weight[0] = M_PI - a;
    color[1]  = clr_2;
    weight[1] = M_PI_2 - weight[0]; 
  }

  if ((a>M_PI) && (a<=1.5*M_PI)){
    color[0]  = clr_2;
    weight[0] = 1.5*M_PI - a;
    color[1]  = clr_3;
    weight[1] = M_PI_2 - weight[0]; 
  }

  if ((a>1.5*M_PI) && (a<=2.0*M_PI)){
    color[0]  = clr_3;
    weight[0] = 2*M_PI - a;
    color[1]  = clr_0;
    weight[1] = M_PI_2 - weight[0]; 
  }

  mean(color.begin(),color.end(),weight.begin(),weight.end(),result);
}


/*

  ##############################
  # class RGB_converter_to_HSV #
  ##############################

*/

template<typename Real>
void
RGB_converter_to_HSV<Real>::
from_RGB_to_cartesian(const double r,
                      const double g,
                      const double b,
                      component_type* const x,
                      component_type* const y,
                      component_type* const z){
  
  double h,s,v;
  RGB_to_HSV(r,g,b,&h,&s,&v);

#ifndef M_PI
  const double theta = 2 * 3.14159265358979323846 * h;
#else
  const double theta = 2 * M_PI * h;
#endif

  const float f = s*v;

  *x = static_cast<component_type>(f * cos(theta));
  *y = static_cast<component_type>(f * sin(theta));
  *z = static_cast<component_type>(v);
}


template<typename Real>
void
RGB_converter_to_HSV<Real>::
from_cartesian_to_RGB(const component_type x,
                      const component_type y,
                      const component_type z,
                      double* const r,
                      double* const g,
                      double* const b){

  const double X = static_cast<double>(x);
  const double Y = static_cast<double>(y);
  const double f = sqrt(X*X + Y*Y);
  const double v = static_cast<double>(z);
  const double s = (v>0.0)?(f/v):(0.0);

#ifndef M_PI
  double theta = atan2(Y,X) / (2 * 3.14159265358979323846);
#else
  double theta = atan2(Y,X) / (2 * M_PI);
#endif

  const double h = (theta>0) ? theta : (theta + 1);
  
  HSV_to_RGB(h,s,v,r,g,b);
}


template<typename Real>
void
RGB_converter_to_HSV<Real>::
RGB_to_HSV(const double r,  const double g,  const double b,
           double* const h, double* const s, double* const v){

  const double min   = std::min(r,std::min(g,b));
  const double max   = std::max(r,std::max(g,b));
  const double delta = max - min;
  
  *v = max; 
    
  if (max!=0)
    *s = delta / max;    
  else {
    *s = 0;
    *h = 0;
    return;
  }
    
  if (delta==0) { 
    *s = 0;
    *h = 0;
    return;
  }

  if(r==max) {
      *h = ( g - b ) / delta;        // between yellow & magenta
  }
  else if( g == max ) {
    *h = 2 + ( b - r ) / delta;    // between cyan & yellow
  }
  else {
    *h = 4 + ( r - g ) / delta;    // between magenta & cyan
  }
    
  *h /= 6;    
  if(*h<0){
    *h += 1;
  }

}
  
template<typename Real>
void
RGB_converter_to_HSV<Real>::
HSV_to_RGB(const double h,const double s,const double v,
           double* const r,double* const g,double* const b){
  
  if(s==0) { // achromatic (grey)
    *r = v;
    *g = v;
    *b = v;
    return;
  }
  
  const int i    = static_cast<int>(floor(6*h));
  const double f = 6*h - i;     
  const double p = v * ( 1 - s );
  const double q = v * ( 1 - s * f );
  const double t = v * ( 1 - s * ( 1 - f ) );
 
  switch(i) {
  case 6:
  case 0:
    *r = v;
    *g = t;
    *b = p;
    break;
    
  case 1:
    *r = q;
    *g = v;
    *b = p;
    break;
    
  case 2:
    *r = p;
    *g = v;
    *b = t;
    break;
    
  case 3:
    *r = p;
    *g = q;
    *b = v;
    break;
    
  case 4:
    *r = t;
    *g = p;
    *b = v;
    break;
    
  default:  
    *r = v;
    *g = p;
    *b = q;
    break;
  }   
}



/*

  ##############################
  # class RGB_converter_to_RGB #
  ##############################

*/


template<typename Real>
void
RGB_converter_to_RGB<Real>::
from_RGB_to_cartesian(const double r,
                      const double g,
                      const double b,
                      component_type* const x,
                      component_type* const y,
                      component_type* const z){
  
  *x = static_cast<component_type>(r);
  *y = static_cast<component_type>(g);
  *z = static_cast<component_type>(b);
}

template<typename Real>
void
RGB_converter_to_RGB<Real>::
from_cartesian_to_RGB(const component_type x,
                      const component_type y,
                      const component_type z,
                      double* const r,
                      double* const g,
                      double* const b){
  
  *r = static_cast<double>(x);
  *g = static_cast<double>(y);
  *b = static_cast<double>(z);
}



/*

  ##############################
  # class RGB_converter_to_XYZ #
  ##############################

*/

template<typename Real>
void
RGB_converter_to_XYZ<Real>::
from_RGB_to_cartesian(const double r,
                      const double g,
                      const double b,
                      component_type* const x,
                      component_type* const y,
                      component_type* const z){
  
  *x = static_cast<component_type>(r*0.412453 + g*0.357580 + b*0.180423);
  *y = static_cast<component_type>(r*0.212671 + g*0.715160 + b*0.072169);
  *z = static_cast<component_type>(r*0.019334 + g*0.119193 + b*0.950227);


}


template<typename Real>
void
RGB_converter_to_XYZ<Real>::
from_cartesian_to_RGB(const component_type x,
                      const component_type y,
                      const component_type z,
                      double* const r,
                      double* const g,
                      double* const b){

  

  const double R =
    static_cast<double>(x*3.240479 + y*-1.537150 + z*-0.498535);
  const double G =
    static_cast<double>(x*-0.969256 + y*1.875992 + z*0.041556);
  const double B =
    static_cast<double>(x*0.055648 + y*-0.204043 + z*1.057311);

  if (R<0.0)      { *r = 0.0; }
  else if (R>1.0) { *r = 1.0; }
  else            { *r = R; }

  if (G<0.0)      { *g = 0.0; }
  else if (G>1.0) { *g = 1.0; }
  else            { *g = G; }

  if (B<0.0)      { *b = 0.0; }
  else if (B>1.0) { *b = 1.0; }
  else            { *b = B; }
}



/*

  ##############################
  # class RGB_converter_to_Luv #
  ##############################
*/

template<typename Real>
const double RGB_converter_to_Luv<Real>::Xn = 0.95045592705;
template<typename Real>
const double RGB_converter_to_Luv<Real>::Yn = 1.0;
template<typename Real>
const double RGB_converter_to_Luv<Real>::Zn = 1.0890577508;

template<typename Real>
double
RGB_converter_to_Luv<Real>::L2Y(const double l){

  return (l>7.99625) ?
    Yn*pow((l+16)/116,3.0) :
    Yn*l/903.3;
}

template<typename Real>
double
RGB_converter_to_Luv<Real>::Y2L(const double y){

  const double ratio = y/Yn;

  return (ratio>0.008856) ?
    116 * pow(ratio,1.0/3.0) - 16 :
    903.3 * ratio;
}


template<typename Real>
void
RGB_converter_to_Luv<Real>::
from_RGB_to_cartesian(const double r,
                      const double g,
                      const double b,
                      component_type* const x,
                      component_type* const y,
                      component_type* const z){
  
  const double X = r*0.412453 + g*0.357580 + b*0.180423;
  const double Y = r*0.212671 + g*0.715160 + b*0.072169;
  const double Z = r*0.019334 + g*0.119193 + b*0.950227;


  *x = static_cast<component_type>(Y2L(Y / Yn));
  const double tmp = Xn + 15 * Yn + 3 * Zn;
  const double unp = 4 * Xn / tmp;
  const double vnp = 9 * Yn / tmp;
  const double tmp2 = X + 15 * Y + 3 * Z;
  const double up = (X==0)? 0 : 4 * X / tmp2;
  const double vp = (Y==0)? 0 : 9 * Y / tmp2;
  *y = static_cast<component_type>(13 * (*x) * (up - unp));
  *z = static_cast<component_type>(13 * (*x) * (vp - vnp));
}


template<typename Real>
void
RGB_converter_to_Luv<Real>::
from_cartesian_to_RGB(const component_type x,
                      const component_type y,
                      const component_type z,
                      double* const r,
                      double* const g,
                      double* const b){

  const double Y = L2Y(static_cast<double>(x));
  const double tmp = Xn + 15 * Yn + 3 * Zn;
  const double unp = 4 * Xn / tmp;
  const double vnp = 9 * Yn / tmp;
  const double Q = ((y==0) ? 0
    : static_cast<double>(y) / (13 * static_cast<double>(x))) + unp;
  const double R = ((z==0) ? 0
    : static_cast<double>(z) / (13 * static_cast<double>(x))) + vnp;
  const double A = 3 * Y * (5 * R - 3);
  const double Z = ((Q - 4) * A - 15 * Q * R * Y) / (12 * R);
  const double X = -(A / R + 3 * Z);


  const double red   = X*3.240479  + Y*-1.537150 + Z*-0.498535;
  const double green = X*-0.969256 + Y*1.875992  + Z*0.041556;
  const double blue  = X*0.055648  + Y*-0.204043 + Z*1.057311;

  if (red<0.0)      { *r = 0.0; }
  else if (red>1.0) { *r = 1.0; }
  else              { *r = red; }

  if (green<0.0)      { *g = 0.0; }
  else if (green>1.0) { *g = 1.0; }
  else                { *g = green; }

  if (blue<0.0)      { *b = 0.0; }
  else if (blue>1.0) { *b = 1.0; }
  else               { *b = blue; }
}


/*

  #####################################
  # class RGB_converter_to_grey_level #
  #####################################
  
*/



template<typename Real>
void
RGB_converter_to_grey_level<Real>::
from_RGB_to_cartesian(const double r,
                      const double g,
                      const double b,
                      component_type* const x,
                      component_type* const y,
                      component_type* const z){

  *x = static_cast<component_type>((r+g+b)/3.0);
  *y = static_cast<component_type>(0);
  *z = static_cast<component_type>(0);
}


template<typename Real>
void
RGB_converter_to_grey_level<Real>::
from_cartesian_to_RGB(const component_type x,
                      const component_type y,
                      const component_type z,
                      double* const r,
                      double* const g,
                      double* const b){
  
  *r = static_cast<double>(x);
  *g = static_cast<double>(x);
  *b = static_cast<double>(x);
}



/*

  #########################
  # Fonctions hors classe #
  #########################

 */

template<typename RGB_converter>
std::ostream& operator<<(std::ostream& s,
                         const Basic_color<RGB_converter>& c){

  s<<c.x<<"\t"<<c.y<<"\t"<<c.z<<"\t";
  
  return s;
}

#undef CONST_STATIC

#endif

---