# I do not do the linear approximation for $u/$un <= 0.008859
sub RGB2Lab {
  my ($red, $green, $blue) = @_;

  my ($x, $y, $z)    = _RGB2XYZitu($red, $green, $blue);
  my ($xn, $yn, $zn) = _RGB2XYZitu(1,1,1);
  my $fx = _cuberoot($x/$xn);
  my $fy = _cuberoot($y/$yn);  
  my $fz = _cuberoot($z/$zn);
  
  return ( 116*$fy-16, 500*($fx - $fy), 200*($fy - $fz) );
 
}

sub Lab2RGB {
  my ($L, $a, $b) = @_;
  my ($xn, $yn, $zn) = _RGB2XYZitu(1,1,1); 
  
  my $fL = _cube(($L+16)/116);
  my $fa = _cube($a/500);
  my $fb = _cube($b/200);
  
  my $y = _cube( ($L+16)/116) * $yn;
  my $x = _cube( ($L+16)/116 + $a/500 ) * $xn;
  my $z = _cube( ($L+16)/116 - $b/200 ) * $zn;
   
  return _XYZitu2RGB($x, $y, $z);

}

sub _RGB2XYZitu {
  my ($r, $g, $b) = @_;
  return ( 
    0.431*$r + 0.342*$g + 0.178*$b,
    0.222*$r + 0.707*$g + 0.071*$b,
    0.020*$r + 0.130*$g + 0.939*$b 
  );
}

sub _XYZitu2RGB {
  my ($x, $y, $z) = @_;
  
  return  map { $_ > 1 ? 1 : $_ } (
    3.063*$x - 1.393*$y - 0.476*$z,
   -0.969*$x + 1.876*$y + 0.042*$z,
    0.068*$x - 0.229*$y + 1.069*$z 
  );

}

sub _cuberoot {
  my $x = shift;  
  return 0 if $x == 0;
  my $sign = ($x < 0) ? -1 : 1;
  $x *= $sign;
  return $sign * exp( log($x)/3.0 );  
}

sub _cube {
  my $x = shift;
  return 0 if $x == 0;
  my $sign = ($x < 0) ? -1 : 1;
  $x *= $sign;
  return $sign * exp( 3 * log($x) );
}