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

#include "vtkArrayDispatch.h"
#include "vtkArrayDispatchDataSetArrayList.h"
#include "vtkDataArrayRange.h"
#include "vtkImageData.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPointSet.h"
#include "vtkPoints.h"
#include "vtkSMPTools.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkUnsignedCharArray.h"

VTK_ABI_NAMESPACE_BEGIN
vtkStandardNewMacro(vtkPointOccupancyFilter);

//------------------------------------------------------------------------------
// Helper classes to support efficient computing, and threaded execution.
namespace
{

//------------------------------------------------------------------------------
// The threaded core of the algorithm. Operator() processes templated points.
template <typename TArray>
struct ComputeOccupancyFunctor
{
  TArray* Points;
  double hX, hY, hZ; // internal data members for performance
  double fX, fY, fZ, bX, bY, bZ;
  vtkIdType xD, yD, zD, xyD;
  unsigned char OccupiedValue;
  unsigned char* Occupancy;

  ComputeOccupancyFunctor(TArray* pts, int dims[3], double origin[3], double spacing[3],
    unsigned char empty, unsigned char occupied, unsigned char* occ)
    : Points(pts)
    , OccupiedValue(occupied)
    , Occupancy(occ)
  {
    std::fill_n(this->Occupancy, dims[0] * dims[1] * dims[2], empty);
    for (int i = 0; i < 3; ++i)
    {
      this->hX = spacing[0];
      this->hY = spacing[1];
      this->hZ = spacing[2];
      this->fX = 1.0 / spacing[0];
      this->fY = 1.0 / spacing[1];
      this->fZ = 1.0 / spacing[2];
      this->bX = origin[0] - 0.5 * this->hX;
      this->bY = origin[1] - 0.5 * this->hY;
      this->bZ = origin[2] - 0.5 * this->hZ;
      this->xD = dims[0];
      this->yD = dims[1];
      this->zD = dims[2];
      this->xyD = dims[0] * dims[1];
    }
  }

  void operator()(vtkIdType ptId, vtkIdType endPtId)
  {
    auto x = vtk::DataArrayTupleRange<3>(this->Points, ptId, endPtId).begin();
    unsigned char* o = this->Occupancy;
    unsigned char ov = this->OccupiedValue;
    int i, j, k;

    for (; ptId < endPtId; ++ptId, ++x)
    {

      i = static_cast<int>((((*x)[0] - this->bX) * this->fX));
      j = static_cast<int>((((*x)[1] - this->bY) * this->fY));
      k = static_cast<int>((((*x)[2] - this->bZ) * this->fZ));

      // If not inside image then skip
      if (i < 0 || i >= this->xD || j < 0 || j >= this->yD || k < 0 || k >= this->zD)
      {
        continue;
      }

      o[i + j * this->xD + k * this->xyD] = ov;

    } // over points
  }
}; // ComputeOccupancyFunctor

struct ComputeOccupancyWorker
{
  template <typename TArray>
  void operator()(TArray* inPts, int dims[3], double origin[3], double spacing[3],
    unsigned char empty, unsigned char occupied, unsigned char* occ)
  {
    ComputeOccupancyFunctor<TArray> functor(inPts, dims, origin, spacing, empty, occupied, occ);
    vtkSMPTools::For(0, inPts->GetNumberOfTuples(), functor);
  }
};
} // anonymous namespace

//================= Begin class proper =======================================
//------------------------------------------------------------------------------
vtkPointOccupancyFilter::vtkPointOccupancyFilter()
{
  this->SampleDimensions[0] = 100;
  this->SampleDimensions[1] = 100;
  this->SampleDimensions[2] = 100;

  // All of these zeros mean automatic computation
  this->ModelBounds[0] = 0.0;
  this->ModelBounds[1] = 0.0;
  this->ModelBounds[2] = 0.0;
  this->ModelBounds[3] = 0.0;
  this->ModelBounds[4] = 0.0;
  this->ModelBounds[5] = 0.0;

  this->Origin[0] = this->Origin[1] = this->Origin[2] = 0.0;
  this->Spacing[0] = this->Spacing[1] = this->Spacing[2] = 1.0;

  this->EmptyValue = 0;
  this->OccupiedValue = 1;
}

//------------------------------------------------------------------------------
vtkPointOccupancyFilter::~vtkPointOccupancyFilter() = default;

//------------------------------------------------------------------------------
int vtkPointOccupancyFilter::FillInputPortInformation(int, vtkInformation* info)
{
  info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPointSet");
  return 1;
}

//------------------------------------------------------------------------------
int vtkPointOccupancyFilter::RequestInformation(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* outputVector)
{
  // get the info objects
  vtkInformation* outInfo = outputVector->GetInformationObject(0);

  int i;
  double ar[3], origin[3];

  outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), 0, this->SampleDimensions[0] - 1,
    0, this->SampleDimensions[1] - 1, 0, this->SampleDimensions[2] - 1);

  for (i = 0; i < 3; i++)
  {
    origin[i] = this->ModelBounds[2 * i];
    if (this->SampleDimensions[i] <= 1)
    {
      ar[i] = 1;
    }
    else
    {
      ar[i] =
        (this->ModelBounds[2 * i + 1] - this->ModelBounds[2 * i]) / (this->SampleDimensions[i] - 1);
    }
  }
  outInfo->Set(vtkDataObject::ORIGIN(), origin, 3);
  outInfo->Set(vtkDataObject::SPACING(), ar, 3);

  vtkDataObject::SetPointDataActiveScalarInfo(outInfo, VTK_UNSIGNED_CHAR, 1);

  return 1;
}

//------------------------------------------------------------------------------
// Compute the size of the sample bounding box automatically from the
// input data.
void vtkPointOccupancyFilter::ComputeModelBounds(
  vtkDataSet* input, vtkImageData* output, vtkInformation* outInfo)
{
  int i;

  // compute model bounds if not set previously
  if (this->ModelBounds[0] >= this->ModelBounds[1] ||
    this->ModelBounds[2] >= this->ModelBounds[3] || this->ModelBounds[4] >= this->ModelBounds[5])
  {
    input->GetBounds(this->ModelBounds);
  }

  // Set volume origin and data spacing
  outInfo->Set(
    vtkDataObject::ORIGIN(), this->ModelBounds[0], this->ModelBounds[2], this->ModelBounds[4]);
  memcpy(this->Origin, outInfo->Get(vtkDataObject::ORIGIN()), sizeof(double) * 3);
  output->SetOrigin(this->Origin);

  for (i = 0; i < 3; i++)
  {
    this->Spacing[i] =
      (this->ModelBounds[2 * i + 1] - this->ModelBounds[2 * i]) / (this->SampleDimensions[i] - 1);
    if (this->Spacing[i] <= 0.0)
    {
      this->Spacing[i] = 1.0;
    }
  }
  outInfo->Set(vtkDataObject::SPACING(), this->Spacing, 3);
  output->SetSpacing(this->Spacing);
}

//------------------------------------------------------------------------------
// Set the dimensions of the sampling volume
void vtkPointOccupancyFilter::SetSampleDimensions(int i, int j, int k)
{
  int dim[3];

  dim[0] = i;
  dim[1] = j;
  dim[2] = k;

  this->SetSampleDimensions(dim);
}

//------------------------------------------------------------------------------
void vtkPointOccupancyFilter::SetSampleDimensions(int dim[3])
{
  int dataDim, i;

  vtkDebugMacro(<< " setting SampleDimensions to (" << dim[0] << "," << dim[1] << "," << dim[2]
                << ")");

  if (dim[0] != this->SampleDimensions[0] || dim[1] != this->SampleDimensions[1] ||
    dim[2] != this->SampleDimensions[2])
  {
    if (dim[0] < 1 || dim[1] < 1 || dim[2] < 1)
    {
      vtkErrorMacro(<< "Bad Sample Dimensions, retaining previous values");
      return;
    }

    for (dataDim = 0, i = 0; i < 3; i++)
    {
      if (dim[i] > 1)
      {
        dataDim++;
      }
    }

    if (dataDim < 3)
    {
      vtkErrorMacro(<< "Sample dimensions must define a volume!");
      return;
    }

    for (i = 0; i < 3; i++)
    {
      this->SampleDimensions[i] = dim[i];
    }

    this->Modified();
  }
}

//------------------------------------------------------------------------------
// Produce the output data
int vtkPointOccupancyFilter::RequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  // get the info objects
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation* outInfo = outputVector->GetInformationObject(0);

  // get the input and output
  vtkPointSet* input = vtkPointSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkImageData* output = vtkImageData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));

  // Check the input
  if (!input || !output)
  {
    return 1;
  }
  vtkIdType numPts = input->GetNumberOfPoints();
  if (numPts < 1)
  {
    return 1;
  }

  // Configure the output
  output->SetExtent(outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()));
  output->AllocateScalars(outInfo);
  int* extent = this->GetExecutive()->GetOutputInformation(0)->Get(
    vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT());

  output->SetDimensions(this->GetSampleDimensions());
  this->ComputeModelBounds(input, output, outInfo);

  // Make sure points are available
  vtkIdType npts = input->GetNumberOfPoints();
  if (npts == 0)
  {
    vtkWarningMacro(<< "No POINTS input!!");
    return 1;
  }

  // Grab the occupancy image and process it.
  output->AllocateScalars(outInfo);
  vtkDataArray* occ = output->GetPointData()->GetScalars();
  unsigned char* o = static_cast<unsigned char*>(output->GetArrayPointerForExtent(occ, extent));

  int dims[3];
  double origin[3], spacing[3];
  output->GetDimensions(dims);
  output->GetOrigin(origin);
  output->GetSpacing(spacing);
  unsigned char ev = this->EmptyValue;
  unsigned char ov = this->OccupiedValue;

  ComputeOccupancyWorker worker;
  if (!vtkArrayDispatch::DispatchByArray<vtkArrayDispatch::PointArrays>::Execute(
        input->GetPoints()->GetData(), worker, dims, origin, spacing, ev, ov, o))
  {
    worker(input->GetPoints()->GetData(), dims, origin, spacing, ev, ov, o);
  }

  return 1;
}

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

  os << indent << "Sample Dimensions: (" << this->SampleDimensions[0] << ", "
     << this->SampleDimensions[1] << ", " << this->SampleDimensions[2] << ")\n";

  os << indent << "ModelBounds: \n";
  os << indent << "  Xmin,Xmax: (" << this->ModelBounds[0] << ", " << this->ModelBounds[1] << ")\n";
  os << indent << "  Ymin,Ymax: (" << this->ModelBounds[2] << ", " << this->ModelBounds[3] << ")\n";
  os << indent << "  Zmin,Zmax: (" << this->ModelBounds[4] << ", " << this->ModelBounds[5] << ")\n";

  os << indent << "Empty Value: " << this->EmptyValue << "\n";
  os << indent << "Occupied Value: " << this->OccupiedValue << "\n";
}
VTK_ABI_NAMESPACE_END