/*=========================================================================

  Program:   Visualization Toolkit

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
#include "vtkLICNoiseHelper.h"

#include "vtkBase64Utilities.h"
#include "vtkGenericDataObjectReader.h"
#include "vtkImageData.h"
#include "vtkNoise200x200.h"

#include <cassert>

namespace
{

/**
integer log base 2
*/
int ilog2(unsigned int n)
{
  if (n == 0)
  {
    return -1;
  }
  unsigned int r = 0;
  while ((n >>= 1))
  {
    r += 1;
  }
  return r;
}

} // end anonymous namespace


//-----------------------------------------------------------------------------
void vtkLICRandomNoise2D::GetValidDimensionAndGrainSize(int type, int &sideLen, int &grainSize)
{
  // perlin noise both side len and grain size need to be powers of 2
  if (type == PERLIN)
  {
    sideLen = 1 << ilog2(sideLen);
    grainSize = 1 << ilog2(grainSize);
  }

  // grains can't be larger than the patch
  if (sideLen < grainSize)
  {
    sideLen = grainSize;
  }

  // generate noise with agiven grainSize size on the patch
  if (sideLen % grainSize)
  {
    // grainSize is not an even divsior of sideLen, adjust sideLen to
    // next larger even divisor
    sideLen = grainSize * (sideLen/grainSize + 1);
  }
}

//-----------------------------------------------------------------------------
int vtkLICRandomNoise2D::ShouldGenerateValue(double prob)
{
  if (this->ProbGen.GetRandomNumber() > (1.0 - prob))
  {
    return 1;
  }
  return 0;
}

//-----------------------------------------------------------------------------
float *vtkLICRandomNoise2D::Generate(
      int type,
      int &sideLen,
      int &grainSize,
      float minNoiseVal,
      float maxNoiseVal,
      int nLevels,
      double impulseProb,
      float impulseBgNoiseVal,
      int seed)
{
  this->GetValidDimensionAndGrainSize(type, sideLen, grainSize);

  switch (type)
  {
    case GAUSSIAN:
      return this->GenerateGaussian(
            sideLen,
            grainSize,
            minNoiseVal,
            maxNoiseVal,
            nLevels,
            impulseProb,
            impulseBgNoiseVal,
            seed);

    case UNIFORM:
      return this->GenerateUniform(
            sideLen,
            grainSize,
            minNoiseVal,
            maxNoiseVal,
            nLevels,
            impulseProb,
            impulseBgNoiseVal,
            seed);

    case PERLIN:
      return this->GeneratePerlin(
            sideLen,
            grainSize,
            minNoiseVal,
            maxNoiseVal,
            nLevels,
            impulseProb,
            impulseBgNoiseVal,
            seed);
  }
  return nullptr;
}

//-----------------------------------------------------------------------------
float *vtkLICRandomNoise2D::GenerateUniform(
      int sideLen,
      int grainSize,
      float minNoiseVal,
      float maxNoiseVal,
      int nLevels,
      double impulseProb,
      float impulseBgNoiseVal,
      int seed)
{
  // generate a patch of single pixel random values
  // with a uniform distribution and fixed number of levels
  nLevels = nLevels < 1 ? 1 : nLevels;
  int maxLevel = nLevels-1;
  float delta = maxLevel != 0 ? 1.0f/maxLevel : 0.0f;
  minNoiseVal = minNoiseVal < 0.0f ? 0.0f : minNoiseVal;
  maxNoiseVal = maxNoiseVal > 1.0f ? 1.0f : maxNoiseVal;
  float noiseRange = maxNoiseVal - minNoiseVal;
  impulseProb = impulseProb < 0.0 ? 0.0 : impulseProb;
  impulseProb = impulseProb > 1.0 ? 1.0 : impulseProb;
  impulseBgNoiseVal = impulseBgNoiseVal < 0.0f ? 0.0f : impulseBgNoiseVal;
  impulseBgNoiseVal = impulseBgNoiseVal > 1.0f ? 1.0f : impulseBgNoiseVal;
  this->ValueGen.SetSeed(seed);
  this->ProbGen.SetSeed(seed);
  const int sdim = sideLen/grainSize;
  const int sdim2 = sdim*sdim;
  float *rvals=(float*)malloc(sdim2*sizeof(float));
  for (int i=0; i<sdim2; ++i)
  {
    rvals[i] = impulseBgNoiseVal;
  }
  for (int j=0; j<sdim; ++j)
  {
     for (int i=0; i<sdim; ++i)
     {
       int idx=j*sdim+i;

       if ((impulseProb == 1.0) || this->ShouldGenerateValue(impulseProb))
       {
         int l = static_cast<int>(this->ValueGen.GetRandomNumber()*nLevels);
         l = l > maxLevel ? maxLevel : l; // needed for 1.0
         rvals[idx] = nLevels == 1 ? maxNoiseVal : minNoiseVal + (l*delta) * noiseRange;
       }
     }
  }

  // map single pixel random values onto a patch of values of
  // the requested grain size
  const int ncomp = 2;
  const int dim2 = sideLen*sideLen;
  const int ntup = ncomp*dim2;
  float *noise = (float*)malloc(ntup*sizeof(float));
  for (int j=0; j<sideLen; ++j)
  {
     for (int i=0; i<sideLen; ++i)
     {
       int idx=ncomp*(j*sideLen+i);

       int ii = i/grainSize;
       int jj = j/grainSize;
       int iidx = jj*sdim+ii;

       noise[idx] = rvals[iidx];
       noise[idx+1] = 1.0f; // alpha
     }
  }
  free(rvals);

  return noise;
}

//-----------------------------------------------------------------------------
float *vtkLICRandomNoise2D::GenerateGaussian(
      int sideLen,
      int grainSize,
      float minNoiseVal,
      float maxNoiseVal,
      int nLevels,
      double impulseProb,
      float impulseBgNoiseVal,
      int seed)
{
  // the distribution becomes Gaussian as N goes to infinity
  const int N = 2048;

  // generate a patch of single pixel random values
  // with a gaussian distribution
  impulseProb = impulseProb < 0.0 ? 0.0 : impulseProb;
  impulseProb = impulseProb > 1.0 ? 1.0 : impulseProb;
  impulseBgNoiseVal = impulseBgNoiseVal < 0.0f ? 0.0f : impulseBgNoiseVal;
  impulseBgNoiseVal = impulseBgNoiseVal > 1.0f ? 1.0f : impulseBgNoiseVal;
  this->ValueGen.SetSeed(seed);
  this->ProbGen.SetSeed(seed);
  const int sdim = sideLen/grainSize;
  const int sdim2 = sdim*sdim;
  float *rvals = (float*)malloc(sdim2*sizeof(float));
  for (int i=0; i<sdim2; ++i)
  {
    rvals[i] = 0.0f;
  }
  for (int j=0; j<sdim; ++j)
  {
    for (int i=0; i<sdim; ++i)
    {
      int idx = j*sdim+i;

      if ((impulseProb == 1.0) || this->ShouldGenerateValue(impulseProb))
      {
        double val = 0.0;
        for (int q=0; q<N; ++q)
        {
          val += this->ValueGen.GetRandomNumber();
        }
        rvals[idx] = static_cast<float>(val);
      }
    }
  }

  // normalize noise field from eps to nLevels onto 0 to 1
  // and restrict to the requested number of levels
  // min/max
  float minVal = static_cast<float>(N+1);
  float maxVal = 0.0f;
  for (int i=0; i<sdim2; ++i)
  {
    // for impulseProb < 1 background is 0 but pixels that are touched
    // have a much larger value, after normalization the gaussian
    // distribution is compressed and localized near 1. We can fix this
    // by ignoring zero values.
    minVal = impulseProb == 1.0 ?
            (rvals[i] < minVal ? rvals[i] : minVal) :
            (rvals[i] < minVal && rvals[i] > 0.0f ? rvals[i] : minVal);

    maxVal = rvals[i]>maxVal ? rvals[i] : maxVal;
  }
  float maxMinDiff = maxVal-minVal;
  // because we ignore zero when impulseProb<1 we have to be careful
  // here so that we can support one noise level.
  minVal = maxMinDiff == 0.0f ? 0.0f : minVal;
  maxMinDiff = maxMinDiff == 0.0f ? (maxVal == 0.0f ? 1.0f : maxVal) : maxMinDiff;

  nLevels = nLevels < 1 ? 1 : nLevels;
  int maxLevel = nLevels-1;
  float delta = maxLevel != 0 ? 1.0f/maxLevel : 0.0f;
  minNoiseVal = minNoiseVal < 0.0f ? 0.0f : minNoiseVal;
  maxNoiseVal = maxNoiseVal > 1.0f ? 1.0f : maxNoiseVal;
  float noiseRange = maxNoiseVal - minNoiseVal;
  for (int i=0; i<sdim2; ++i)
  {
    // normalize
    float val = rvals[i] < minVal ? rvals[i] : (rvals[i] - minVal)/maxMinDiff;
    // restrict
    int l = static_cast<int>(val*nLevels);
    l = l > maxLevel ? maxLevel : l;
    rvals[i]
       = rvals[i] < minVal ? impulseBgNoiseVal
       : nLevels == 1 ? maxNoiseVal : minNoiseVal + (l*delta) * noiseRange;
  }

  // map single pixel random values onto a patch of values of
  // the requested grain size
  const int ncomp = 2;
  const int dim2 = sideLen*sideLen;
  const int ntup = ncomp*dim2;
  float *noise = (float*)malloc(ntup*sizeof(float));
  for (int j=0; j<sideLen; ++j)
  {
     for (int i=0; i<sideLen; ++i)
     {
       int idx = ncomp*(j*sideLen+i);

       int ii = i/grainSize;
       int jj = j/grainSize;
       int iidx = jj*sdim+ii;

       noise[idx] = rvals[iidx];
       noise[idx+1] = 1.0; // alpha
     }
  }
  free(rvals);

  return noise;
}

//-----------------------------------------------------------------------------
float *vtkLICRandomNoise2D::GeneratePerlin(
      int sideLen,
      int grainSize,
      float minNoiseVal,
      float maxNoiseVal,
      int nLevels,
      double impulseProb,
      float impulseBgNoiseVal,
      int seed)
{
  // note: requires power of 2 sideLen, and sideLen > grainSize
  const int ncomp = 2;
  const int dim2 = sideLen*sideLen;
  const int ntup = ncomp*dim2;
  float *noise = static_cast<float*>(malloc(ntup*sizeof(float)));
  for (int i=0; i<ntup; i+=2)
  {
    noise[i  ] = 0.0f;
    noise[i+1] = 1.0f; // alpha channel
  }

  impulseProb = impulseProb < 0.0 ? 0.0 : impulseProb;
  impulseProb = impulseProb > 1.0 ? 1.0 : impulseProb;
  impulseBgNoiseVal = impulseBgNoiseVal < 0.0f ? 0.0f : impulseBgNoiseVal;
  impulseBgNoiseVal = impulseBgNoiseVal > 1.0f ? 1.0f : impulseBgNoiseVal;
  minNoiseVal = minNoiseVal < 0.0f ? 0.0f : minNoiseVal;
  maxNoiseVal = maxNoiseVal > 1.0f ? 1.0f : maxNoiseVal;

  //int nIter = ilog2(static_cast<unsigned int>(sideLen-1<nLevels ? sideLen-1 : nLevels));
  int nIter = ilog2(static_cast<unsigned int>(grainSize));
  for (int w=0; w<nIter; ++w)
  {
    // reduce range with grain size
    float levelNoiseMin = 0.0f;
    float levelNoiseMax = 0.1f + 0.9f/static_cast<float>(1<<(nIter-1-w));
    //float levelNoiseMax = 1.0f - levelNoiseMin;
    // generate a level of noise
    int levelGrainSize = 1<<w;
    float *levelNoise = GenerateGaussian(
          sideLen,
          levelGrainSize,
          levelNoiseMin,
          levelNoiseMax,
          nLevels,
          impulseProb,
          impulseBgNoiseVal,
          seed);
    /*// smooth
    int nsp = w;
    for (int k=0; k<nsp; ++k)
      {
      for (int j=0; j<sideLen; ++j)
        {
         for (int i=0; i<sideLen; ++i)
           {
           float K[9] = {
             0.0191724, 0.100120, 0.0191724,
             0.1001200, 0.522831, 0.1001200,
             0.0191724, 0.100120, 0.0191724
             };
           float val=0.0;
           for (int q=0; q<3; ++q)
             {
             for (int p=0; p<3; ++p)
               {
               int ii = i+p-1;
               ii = ii < 0 ? i : ii;
               ii = ii >= sideLen ? i : ii;
               int jj = j+q-1;
               jj = jj < 0 ? j : jj;
               jj = jj >= sideLen ? j : jj;
               int idx = 2*(sideLen*jj+ii);
               val += levelNoise[idx]*K[q*3+p];
               }
             }
           levelNoise[2*(sideLen*j+i)] = val;
           }
        }
      }*/
    // accumulate
    for (int i=0; i<ntup; i+=2)
    {
      noise[i] += levelNoise[i];
    }
    free(levelNoise);
  }
  // normalize
  float minVal = static_cast<float>(nIter+1);
  float maxVal = 0.0f;
  for (int i=0; i<ntup; i+=2)
  {
    float val = noise[i];
    minVal = val<minVal ? val : minVal;
    maxVal = val>maxVal ? val : maxVal;
  }
  float maxMinDiff = maxVal - minVal;
  if ( maxMinDiff <= 0.0f )
  {
    maxMinDiff = 1.0f;
    minVal = 0.0f;
  }
  for (int i=0; i<ntup; i+=2)
  {
    noise[i] = (noise[i] - minVal) / maxMinDiff;
  }
  return noise;
}

/**
Load a predefiined texture that has been "pickled" in a string.
This texture is 200x200 pixles, has a Gaussian distribution, and
intensities ranging between 0 and 206. This is the texture that
is used when GenerateNoiseTexture is disabled.
*/
vtkImageData *vtkLICRandomNoise2D::GetNoiseResource()
{
  std::string base64string;
  for (unsigned int cc=0; cc < file_noise200x200_vtk_nb_sections; cc++)
  {
    base64string += reinterpret_cast<const char*>(file_noise200x200_vtk_sections[cc]);
  }

  unsigned char* binaryInput
     = new unsigned char[file_noise200x200_vtk_decoded_length + 10];

  unsigned long binarylength = static_cast<unsigned long>(
    vtkBase64Utilities::DecodeSafely(
        reinterpret_cast<const unsigned char*>(base64string.c_str()),
        base64string.length(),
        binaryInput,
        file_noise200x200_vtk_decoded_length + 10));

  assert("check valid_length"
    && (binarylength == file_noise200x200_vtk_decoded_length));

  vtkGenericDataObjectReader* reader = vtkGenericDataObjectReader::New();
  reader->ReadFromInputStringOn();

  reader->SetBinaryInputString(
        reinterpret_cast<char*>(binaryInput),
        static_cast<int>(binarylength));

  reader->Update();
  vtkImageData* data = vtkImageData::New();
  data->ShallowCopy(reader->GetOutput());

  delete [] binaryInput;
  reader->Delete();
  return data;
}