// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause
#include "vtkEncodedGradientEstimator.h"

#include "vtkGarbageCollector.h"
#include "vtkImageData.h"
#include "vtkMultiThreader.h"
#include "vtkRecursiveSphereDirectionEncoder.h"
#include "vtkTimerLog.h"

#include <cmath>

VTK_ABI_NAMESPACE_BEGIN
vtkCxxSetObjectMacro(vtkEncodedGradientEstimator, InputData, vtkImageData);

// Construct a vtkEncodedGradientEstimator with initial values of nullptr for
// the Input, EncodedNormal, and GradientMagnitude. Also,
// indicate that the IndexTable has not yet been initialized. The
// GradientMagnitudeRange and the GradientMangitudeTable are
// initialized to default values - these will change in the future
// when magnitude of gradient opacities are included
vtkEncodedGradientEstimator::vtkEncodedGradientEstimator()
{
  this->InputData = nullptr;
  this->EncodedNormals = nullptr;
  this->EncodedNormalsSize[0] = 0;
  this->EncodedNormalsSize[1] = 0;
  this->EncodedNormalsSize[2] = 0;
  this->GradientMagnitudes = nullptr;
  this->GradientMagnitudeScale = 1.0;
  this->GradientMagnitudeBias = 0.0;
  this->Threader = vtkMultiThreader::New();
  this->NumberOfThreads = this->Threader->GetNumberOfThreads();
  this->DirectionEncoder = vtkRecursiveSphereDirectionEncoder::New();
  this->ComputeGradientMagnitudes = 1;
  this->CylinderClip = 0;
  this->CircleLimits = nullptr;
  this->CircleLimitsSize = -1;
  this->UseCylinderClip = 0;
  this->LastUpdateTimeInSeconds = -1.0;
  this->LastUpdateTimeInCPUSeconds = -1.0;
  this->ZeroNormalThreshold = 0.0;
  this->ZeroPad = 1;
  this->BoundsClip = 0;
  this->Bounds[0] = this->Bounds[1] = this->Bounds[2] = this->Bounds[3] = this->Bounds[4] =
    this->Bounds[5] = 0;
}

// Destruct a vtkEncodedGradientEstimator - free up any memory used
vtkEncodedGradientEstimator::~vtkEncodedGradientEstimator()
{
  this->SetInputData(nullptr);
  this->Threader->Delete();
  this->Threader = nullptr;

  delete[] this->EncodedNormals;

  delete[] this->GradientMagnitudes;

  if (this->DirectionEncoder)
  {
    this->DirectionEncoder->UnRegister(this);
  }

  delete[] this->CircleLimits;
}

void vtkEncodedGradientEstimator::SetZeroNormalThreshold(float v)
{
  if (this->ZeroNormalThreshold != v)
  {
    if (v < 0.0)
    {
      vtkErrorMacro(<< "The ZeroNormalThreshold must be a value >= 0.0");
      return;
    }

    this->ZeroNormalThreshold = v;
    this->Modified();
  }
}

void vtkEncodedGradientEstimator::SetDirectionEncoder(vtkDirectionEncoder* direnc)
{
  // If we are setting it to its current value, don't do anything
  if (this->DirectionEncoder == direnc)
  {
    return;
  }

  // If we already have a direction encoder, unregister it.
  if (this->DirectionEncoder)
  {
    this->DirectionEncoder->UnRegister(this);
    this->DirectionEncoder = nullptr;
  }

  // If we are passing in a non-nullptr encoder, register it
  if (direnc)
  {
    direnc->Register(this);
  }

  // Actually set the encoder, and consider the object Modified
  this->DirectionEncoder = direnc;
  this->Modified();
}

int vtkEncodedGradientEstimator::GetEncodedNormalIndex(vtkIdType xyzIndex)
{
  this->Update();
  return *(this->EncodedNormals + xyzIndex);
}

int vtkEncodedGradientEstimator::GetEncodedNormalIndex(int xIndex, int yIndex, int zIndex)
{
  vtkIdType ystep, zstep;

  this->Update();

  // Compute steps through the volume in x, y, and z
  ystep = this->InputSize[0];
  zstep = ystep * this->InputSize[1];

  return *(this->EncodedNormals + zIndex * zstep + yIndex * ystep + xIndex);
}

unsigned short* vtkEncodedGradientEstimator::GetEncodedNormals()
{
  this->Update();

  return this->EncodedNormals;
}

unsigned char* vtkEncodedGradientEstimator::GetGradientMagnitudes()
{
  this->Update();

  return this->GradientMagnitudes;
}

void vtkEncodedGradientEstimator::Update()
{
  int scalarInputSize[3];
  double scalarInputAspect[3];
  double startSeconds, endSeconds;
  double startCPUSeconds, endCPUSeconds;

  if (!this->InputData)
  {
    vtkErrorMacro(<< "No input in gradient estimator.");
    return;
  }

  if (this->GetMTime() > this->BuildTime || this->DirectionEncoder->GetMTime() > this->BuildTime ||
    this->InputData->GetMTime() > this->BuildTime || !this->EncodedNormals)
  {

    startSeconds = vtkTimerLog::GetUniversalTime();
    startCPUSeconds = vtkTimerLog::GetCPUTime();

    // Get the dimensions of the data and its aspect ratio
    this->InputData->GetDimensions(scalarInputSize);
    this->InputData->GetSpacing(scalarInputAspect);

    // If we previously have allocated space for the encoded normals,
    // and this space is no longer the right size, delete it
    if (this->EncodedNormalsSize[0] != scalarInputSize[0] ||
      this->EncodedNormalsSize[1] != scalarInputSize[1] ||
      this->EncodedNormalsSize[2] != scalarInputSize[2])
    {
      delete[] this->EncodedNormals;
      this->EncodedNormals = nullptr;

      delete[] this->GradientMagnitudes;
      this->GradientMagnitudes = nullptr;
    }

    // Compute the number of encoded voxels
    vtkIdType encodedSize = scalarInputSize[0];
    encodedSize *= scalarInputSize[1];
    encodedSize *= scalarInputSize[2];

    // Allocate space for the encoded normals if necessary
    if (!this->EncodedNormals)
    {
      this->EncodedNormals = new unsigned short[encodedSize];
      this->EncodedNormalsSize[0] = scalarInputSize[0];
      this->EncodedNormalsSize[1] = scalarInputSize[1];
      this->EncodedNormalsSize[2] = scalarInputSize[2];
    }

    if (!this->GradientMagnitudes && this->ComputeGradientMagnitudes)
    {
      this->GradientMagnitudes = new unsigned char[encodedSize];
    }

    // Copy info that multi threaded function will need into temp variables
    memcpy(this->InputSize, scalarInputSize, 3 * sizeof(int));
    // TODO cleanup when double changes are further along
    this->InputAspect[0] = static_cast<float>(scalarInputAspect[0]);
    this->InputAspect[1] = static_cast<float>(scalarInputAspect[1]);
    this->InputAspect[2] = static_cast<float>(scalarInputAspect[2]);
    // memcpy( this->InputAspect, scalarInputAspect, 3 * sizeof(float) );

    if (this->CylinderClip && (this->InputSize[0] == this->InputSize[1]))
    {
      this->UseCylinderClip = 1;
      this->ComputeCircleLimits(this->InputSize[0]);
    }
    else
    {
      this->UseCylinderClip = 0;
    }
    this->UpdateNormals();

    this->BuildTime.Modified();

    endSeconds = vtkTimerLog::GetUniversalTime();
    endCPUSeconds = vtkTimerLog::GetCPUTime();

    this->LastUpdateTimeInSeconds = static_cast<float>(endSeconds - startSeconds);
    this->LastUpdateTimeInCPUSeconds = static_cast<float>(endCPUSeconds - startCPUSeconds);
  }
}

void vtkEncodedGradientEstimator::ComputeCircleLimits(int size)
{
  int *ptr, y;
  double w, halfsize, length, start, end;

  if (this->CircleLimitsSize != size)
  {
    delete[] this->CircleLimits;
    this->CircleLimits = new int[2 * size];
    this->CircleLimitsSize = size;
  }

  ptr = this->CircleLimits;

  halfsize = (size - 1) / 2.0;

  for (y = 0; y < size; y++)
  {
    w = halfsize - y;
    length = static_cast<int>(sqrt((halfsize * halfsize) - (w * w)) + 0.5);
    start = halfsize - length - 1;
    end = halfsize + length + 1;
    start = (start < 0) ? (0) : (start);
    end = (end > (size - 1)) ? (size - 1) : (end);

    *(ptr++) = static_cast<int>(start);
    *(ptr++) = static_cast<int>(end);
  }
}

// Print the vtkEncodedGradientEstimator
void vtkEncodedGradientEstimator::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  if (this->InputData)
  {
    os << indent << "InputData: (" << this->InputData << ")\n";
  }
  else
  {
    os << indent << "Input: (none)\n";
  }

  if (this->DirectionEncoder)
  {
    os << indent << "DirectionEncoder: (" << this->DirectionEncoder << ")\n";
  }
  else
  {
    os << indent << "DirectionEncoder: (none)\n";
  }

  os << indent << "Build Time: " << this->BuildTime.GetMTime() << endl;

  os << indent << "Gradient Magnitude Scale: " << this->GradientMagnitudeScale << endl;

  os << indent << "Gradient Magnitude Bias: " << this->GradientMagnitudeBias << endl;

  os << indent << "Zero Pad: " << ((this->ZeroPad) ? "On" : "Off") << endl;

  os << indent << "Bounds Clip: " << ((this->BoundsClip) ? "On" : "Off") << endl;

  os << indent << "Bounds: (" << this->Bounds[0] << ", " << this->Bounds[1] << ", "
     << this->Bounds[2] << ", " << this->Bounds[3] << ", " << this->Bounds[4] << ", "
     << this->Bounds[5] << ")\n";

  os << indent << "Zero Normal Threshold: " << this->ZeroNormalThreshold << endl;

  os << indent
     << "Compute Gradient Magnitudes: " << ((this->ComputeGradientMagnitudes) ? "On" : "Off")
     << endl;

  os << indent << "Cylinder Clip: " << ((this->CylinderClip) ? "On" : "Off") << endl;

  os << indent << "Number Of Threads: " << this->NumberOfThreads << endl;

  os << indent << "Last Update Time In Seconds: " << this->LastUpdateTimeInSeconds << endl;

  os << indent << "Last Update Time In CPU Seconds: " << this->LastUpdateTimeInCPUSeconds << endl;

  // I don't want to print out these variables - they are
  // internal and the get methods are included only for access
  // within the threaded function
  // os << indent << "Use Cylinder Clip: "
  //    << this->UseCylinderClip << endl;
  // os << indent << " Input Size: "
  //    << this->InputSize << endl;
  // os << indent << " Input Aspect Clip: "
  //    << this->InputAspect << endl;
}

//------------------------------------------------------------------------------
void vtkEncodedGradientEstimator::ReportReferences(vtkGarbageCollector* collector)
{
  this->Superclass::ReportReferences(collector);
  vtkGarbageCollectorReport(collector, this->InputData, "Input");
}
VTK_ABI_NAMESPACE_END
