// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause

// Description
// This test creates a single sphere inside a volume and fills the sphere
// with a segmented regions id = 1. It then produces sampled points from the
// volume.

#include "vtkActor.h"
#include "vtkFeatureEdges.h"
#include "vtkGeneralizedSurfaceNets3D.h"
#include "vtkIdTypeArray.h"
#include "vtkImageData.h"
#include "vtkIntArray.h"
#include "vtkJitterPoints.h"
#include "vtkLabeledImagePointSampler.h"
#include "vtkMath.h"
#include "vtkNew.h"
#include "vtkOutlineFilter.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkPolyData.h"
#include "vtkPolyDataMapper.h"
#include "vtkPolyDataWriter.h"
#include "vtkProperty.h"
#include "vtkRegressionTestImage.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkRenderer.h"
#include "vtkSurfaceNets3D.h"
#include "vtkTesting.h"
#include "vtkTimerLog.h"
#include "vtkVoronoi3D.h"
#include "vtkWindowedSincPolyDataFilter.h"

namespace
{

// Sphere of radius R centered at (0,0,0)
void LabelVoxels(vtkImageData* volume, double R)
{
  int dims[3];
  volume->GetDimensions(dims);
  vtkIdType numPts = volume->GetNumberOfPoints();
  assert(numPts == (dims[0] * dims[1] * dims[2]));

  vtkDataArray* scalars = volume->GetPointData()->GetScalars();
  int* s = vtkIntArray::FastDownCast(scalars)->GetPointer(0);

  double d2, x[3], R2 = R * R;
  double origin[3] = {0,0,0};
  for (vtkIdType ptId = 0; ptId < numPts; ++ptId)
  {
    volume->GetPoint(ptId, x);
    d2 = vtkMath::Distance2BetweenPoints(x, origin);
    *s++ = (d2 <= R2 ? 1 : 0);
  }
}

} // anonymous

int TestLabeledImagePointSampler(int argc, char* argv[])
{
  // Create a volume
  int DIM = 5;
  double RADIUS = DIM / 2.5;
  vtkNew<vtkImageData> volume;
  volume->SetDimensions(DIM, DIM, DIM);
  volume->SetOrigin(-DIM/2,-DIM/2,-DIM/2);
  //  volume->SetOrigin(-DIM / 2, -DIM / 2, -DIM / 2);
  volume->SetSpacing(1, 1, 1);
  volume->AllocateScalars(VTK_INT, 1);

  // Label the sphere region
  LabelVoxels(volume, RADIUS);

  // Surface net for context
  vtkNew<vtkSurfaceNets3D> snet;
  snet->SetInputData(volume);
  snet->SetLabel(0, 1);

  vtkNew<vtkPolyDataMapper> snetMapper;
  snetMapper->SetInputConnection(snet->GetOutputPort());
  snetMapper->ScalarVisibilityOff();

  vtkNew<vtkActor> snetActor;
  snetActor->SetMapper(snetMapper);
  snetActor->GetProperty()->SetColor(0.5, 0.5, 0.5);
  snetActor->GetProperty()->SetOpacity(0.4);

  // Time processing
  vtkNew<vtkTimerLog> timer;

  // Sampled points
  vtkNew<vtkLabeledImagePointSampler> sampler;
  sampler->SetInputData(volume);
  sampler->SetLabel(0, 1);
  sampler->SetDensityDistributionToExponential();
  sampler->SetDensityDistributionToLinear();
  sampler->SetN(1);
  sampler->SetOutputTypeToLabeledPoints();
  sampler->SetOutputTypeToBackgroundPoints();
  sampler->SetOutputTypeToAllPoints();
  sampler->GenerateVertsOn();
  sampler->RandomizeOff();
  sampler->JitterOff();

  timer->StartTimer();
  sampler->Update();
  timer->StopTimer();
  auto time = timer->GetElapsedTime();
  cout << "Time to sample data: " << time << "\n";

  // Jitter points if desired
  vtkNew<vtkJitterPoints> jitter;
  jitter->SetInputConnection(sampler->GetOutputPort());
  jitter->SetRadius(0.1);
  jitter->RadiusIsAbsoluteOn();
  jitter->SetJitterToUnconstrained();

  timer->StartTimer();
  jitter->Update();
  timer->StopTimer();
  time = timer->GetElapsedTime();
  cout << "Time to jitter points: " << time << "\n";

  vtkNew<vtkIdTypeArray> poi;
  //  poi->InsertNextValue(10);
  poi->InsertNextValue(10);

  // Voronoi-based surface net
  // vtkNew<vtkGeneralizedSurfaceNets3D> vsn;
  // vsn->SetInputConnection(sampler->GetOutputPort());
  // //  vsn->SetInputConnection(jitter->GetOutputPort());
  // vsn->SetLabel(0, 1);
  // vsn->GetLocator()->SetNumberOfPointsPerBucket(2);

  vtkNew<vtkVoronoi3D> vsn;
  vsn->SetInputConnection(sampler->GetOutputPort());
  //  vsn->SetInputConnection(jitter->GetOutputPort());
  //  vsn->SetPointsOfInterest(poi);
  // v->SetMaximumNumberOfHullClips(7);
  vsn->GetLocator()->SetNumberOfPointsPerBucket(2);
  vsn->SetOutputTypeToSurfaceNet();
  vsn->SetOutputTypeToBoundary();
  vsn->ValidateOn();

  timer->StartTimer();
  vsn->Update();
  timer->StopTimer();
  time = timer->GetElapsedTime();
  cout << "Time to surface data: " << time << "\n";
  cout << "\tNumber of threads used: " << vsn->GetNumberOfThreadsUsed() << "\n";

  // Adjacency graph
  vtkNew<vtkVoronoi3D> ag;
  ag->SetInputConnection(sampler->GetOutputPort());
  //  ag->SetInputConnection(jitter->GetOutputPort());
  ag->GetLocator()->SetNumberOfPointsPerBucket(2);
  ag->SetOutputTypeToAdjacencyGraph();

  vtkNew<vtkPolyDataMapper> agMapper;
  agMapper->SetInputConnection(ag->GetOutputPort());

  vtkNew<vtkActor> agActor;
  agActor->SetMapper(agMapper);
  agActor->GetProperty()->SetColor(1,1,1);

  // Smoothing merged points
  vtkNew<vtkWindowedSincPolyDataFilter> smoother;
  smoother->SetInputConnection(vsn->GetOutputPort());
  smoother->SetNumberOfIterations(80);
  smoother->FeatureEdgeSmoothingOff();
  smoother->SetPassBand(0.05);
  // vtkNew<vtkConstrainedSmoothingFilter> smoother;
  // smoother->SetInputConnection(vsn->GetOutputPort());
  // smoother->SetNumberOfIterations(40);
  // smoother->SetConstraintDistance(1);
  // smoother->SetRelaxationFactor(0.5);

  // Edges
  vtkNew<vtkFeatureEdges> feats;
  feats->SetInputConnection(smoother->GetOutputPort());
  feats->ExtractAllEdgeTypesOff();
  feats->BoundaryEdgesOn();
  feats->NonManifoldEdgesOn();

  vtkNew<vtkPolyDataMapper> featsMapper;
  featsMapper->SetInputConnection(feats->GetOutputPort());
  featsMapper->ScalarVisibilityOff();

  vtkNew<vtkActor> featsActor;
  featsActor->SetMapper(featsMapper);
  featsActor->GetProperty()->SetColor(0, 1, 0);
  featsActor->GetProperty()->SetEdgeWidth(2);

  vtkNew<vtkPolyDataMapper> mapper;
  mapper->SetInputConnection(jitter->GetOutputPort());
  mapper->SetInputConnection(vsn->GetOutputPort());
  //  mapper->SetInputConnection(smoother->GetOutputPort());
  mapper->SetScalarRange(0, 64);
  // mapper->SetScalarRange(0,sampler->GetOutput()->GetNumberOfPoints());
  //  mapper->SetScalarRange(0,v->GetOutput()->GetNumberOfCells());
  //  mapper->ScalarVisibilityOff();

  vtkNew<vtkActor> actor;
  actor->SetMapper(mapper);
  actor->GetProperty()->SetColor(1, 1, 1);
  actor->GetProperty()->SetPointSize(2);

  // Report some stats
  sampler->Update();
  cout << "Number of sampled points: " << sampler->GetOutput()->GetNumberOfPoints() << "\n";
  cout << "Number of Voronoi points: " << vsn->GetOutput()->GetNumberOfPoints() << "\n";
  cout << "Number of Voronoi cells: " << vsn->GetOutput()->GetNumberOfCells() << "\n";

  // Bounding box
  vtkNew<vtkOutlineFilter> outline;
  outline->SetInputData(volume);

  vtkNew<vtkPolyDataMapper> outlineMapper;
  outlineMapper->SetInputConnection(outline->GetOutputPort());

  vtkNew<vtkActor> outlineActor;
  outlineActor->SetMapper(outlineMapper);

  vtkNew<vtkRenderer> ren;
  //  ren->AddActor(snetActor);
  ren->AddActor(actor);
  ren->AddActor(agActor);
  ren->AddActor(featsActor);
  //  ren->AddActor(outlineActor);

  vtkNew<vtkRenderWindow> renWin;
  renWin->SetSize(400, 400);
  renWin->AddRenderer(ren);

  vtkNew<vtkRenderWindowInteractor> iren;
  iren->SetRenderWindow(renWin);

  ren->ResetCamera();
  renWin->Render();

  int retVal = vtkRegressionTestImage(renWin);
  if (retVal == vtkRegressionTester::DO_INTERACTOR)
  {
    iren->Start();
  }
  return !retVal;
}
