template<typename MatType, typename CubeType>
CubeType StructuredForests<MatType, CubeType>::
RGB2LUV(CubeType const &InImage)
{
  //assert type is double or float.
  double a, y0, maxi;
  a = std::pow(29.0, 3) / 27.0;
  y0 = 8.0 / a;
  maxi = 1.0 / 270.0;

  arma::vec table(1064);
  for (size_t i = 0; i <= 1024; ++i)
  {
    table(i) = i / 1024.0;

    if (table(i) > y0)
      table(i) = 116 * pow(table(i), 1.0/3.0) - 16.0;
    else
      table(i) = table(i) * a;

    table(i) = table(i) * maxi;
  }
  for(size_t i = 1025; i < table.n_elem; ++i)
  {
    table(i) = table(i - 1);
  }

  MatType rgb2xyz;
  rgb2xyz << 0.430574 << 0.222015 << 0.020183 << arma::endr
          << 0.341550 << 0.706655 << 0.129553 << arma::endr
          << 0.178325 << 0.071330 << 0.939180;

  //see how to calculate this efficiently. numpy.dot does this.
  CubeType xyz(InImage.n_rows, InImage.n_cols, rgb2xyz.n_cols);
  /*
  for(size_t i = 0; i < InImage.n_rows; ++i)
  {
    for(size_t j = 0; j < InImage.n_cols; ++j)
    {
      for(size_t k = 0; k < rgb2xyz.n_cols; ++k)
      {
        double val = 0;
        for(size_t l = 0; l < InImage.n_slices; ++l)
        {
          val = val + InImage(i, j, l) *\
              rgb2xyz(l, k);
        }
        xyz(i, j, k) = val;
      }
    }
  }
  */
  for (size_t i = 0; i < InImage.slice(0).n_elem; ++i)
  {
    double r = InImage.slice(0)(i);
    double g = InImage.slice(1)(i);
    double b = InImage.slice(2)(i);

    xyz.slice(0)(i) = 0.430574 * r + 0.341550 * g + 0.178325 * b;
    xyz.slice(1)(i) = 0.222015 * r + 0.706655 * g + 0.129553 * b;
    xyz.slice(2)(i) = 0.020183 * r + 0.071330 * g + 0.939180 * b;
  }

  std::cout << "printing xyz" << std::endl;
  xyz.print();
  MatType nz(InImage.n_rows, InImage.n_cols);

  nz = 1.0 / ( xyz.slice(0) + (15 * xyz.slice(1) ) +
       (3 * xyz.slice(2) + 1e-35));
  std::cout << arma::size(nz) << std::endl;
  nz.print();
  CubeType OutImage(InImage.n_rows, InImage.n_cols, InImage.n_slices);

  for(size_t j = 0; j < xyz.n_cols; ++j)
  {
    for(size_t i = 0; i < xyz.n_rows; ++i)
    {
      OutImage(i, j, 0) = table( static_cast<size_t>( (1024 * xyz(i, j, 1) ) ) );
    }
  }

  OutImage.slice(1) = OutImage.slice(0) % (13 * 4 * (xyz.slice(0) % nz) \
                        - 13 * 0.197833) + 88 * maxi;
  OutImage.slice(2) = OutImage.slice(0) % (13 * 9 * (xyz.slice(1) % nz) \
                        - 13 * 0.468331) + 134 * maxi;

  return OutImage;
}