VectorFieldNonZeroExtraction

// #include <vtkDataArray.h>
// #include <vtkFloatArray.h>
#include <vtkImageData.h>
#include <vtkImageMagnitude.h>
#include <vtkNew.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkThresholdPoints.h>
#include <vtkXMLImageDataWriter.h>
#include <vtkXMLPolyDataWriter.h>

namespace {
void CreateVectorField(vtkImageData* image);
void WriteImage(vtkImageData* image, const std::string& fileName);
void WriteVectorField(vtkPolyData* vectorField, const std::string& fileName);
} // namespace

int main(int, char*[])
{
  // Create an image
  vtkNew<vtkImageData> image;
  CreateVectorField(image);

  vtkNew<vtkImageMagnitude> magnitudeFilter;
  magnitudeFilter->SetInputData(image);
  magnitudeFilter->Update(); // This filter produces a vtkImageData with an
                             // array named "Magnitude".

  image->GetPointData()->AddArray(
      magnitudeFilter->GetOutput()->GetPointData()->GetScalars());
  image->GetPointData()->SetActiveScalars("Magnitude");

  vtkNew<vtkThresholdPoints> thresholdPoints;
  // thresholdPoints->SetInputConnection(magnitudeFilter->GetOutputPort());
  thresholdPoints->SetInputData(image);
  thresholdPoints->SetUpperThreshold(0.05);
  thresholdPoints->SetThresholdFunction(vtkThresholdPoints::THRESHOLD_UPPER);
  thresholdPoints->Update();

  WriteImage(image, "input.vti");
  WriteVectorField(thresholdPoints->GetOutput(), "output.vtp");

  return EXIT_SUCCESS;
}

namespace {

void CreateVectorField(vtkImageData* image)
{
  // Specify the size of the image data
  image->SetDimensions(50, 50, 1);
  image->AllocateScalars(VTK_FLOAT, 3);

  const int* dims = image->GetDimensions();

  // Zero the image
  for (auto y = 0; y < dims[1]; y++)
  {
    for (auto x = 0; x < dims[0]; x++)
    {
      float* pixel = static_cast<float*>(image->GetScalarPointer(x, y, 0));
      pixel[0] = 0.0;
      pixel[1] = 0.0;
      pixel[2] = 0.0;
    }
  }

  // Set two of the pixels to non zero values
  float* pixel = static_cast<float*>(image->GetScalarPointer(20, 20, 0));
  pixel[0] = -10.0;
  pixel[1] = 5.0;
  pixel[2] = 0.0;

  pixel = static_cast<float*>(image->GetScalarPointer(30, 30, 0));
  pixel[0] = 10.0;
  pixel[1] = 10.0;
  pixel[2] = 0.0;

  image->GetPointData()->SetActiveVectors("ImageScalars");

  image->Modified();
}

void WriteVectorField(vtkPolyData* vectorField, const std::string& fileName)
{
  vtkNew<vtkXMLPolyDataWriter> writer;
  writer->SetFileName(fileName.c_str());
  writer->SetInputData(vectorField);
  writer->Write();
}

void WriteImage(vtkImageData* image, const std::string& fileName)
{
  vtkNew<vtkXMLImageDataWriter> writer;
  writer->SetFileName(fileName.c_str());
  writer->SetInputData(image);
  writer->Write();
}

} // namespace

cmake_minimum_required(VERSION 3.12 FATAL_ERROR)

project(VectorFieldNonZeroExtraction)

find_package(VTK COMPONENTS 
  CommonCore
  CommonDataModel
  FiltersCore
  IOXML
  ImagingMath
)

if (NOT VTK_FOUND)
  message(FATAL_ERROR "VectorFieldNonZeroExtraction: Unable to find the VTK build folder.")
endif()

# Prevent a "command line is too long" failure in Windows.
set(CMAKE_NINJA_FORCE_RESPONSE_FILE "ON" CACHE BOOL "Force Ninja to use response files.")
add_executable(VectorFieldNonZeroExtraction MACOSX_BUNDLE VectorFieldNonZeroExtraction.cxx )
  target_link_libraries(VectorFieldNonZeroExtraction PRIVATE ${VTK_LIBRARIES}
)
# vtk_module_autoinit is needed
vtk_module_autoinit(
  TARGETS VectorFieldNonZeroExtraction
  MODULES ${VTK_LIBRARIES}
)

cd VectorFieldNonZeroExtraction/build

cmake ..

cmake -DVTK_DIR:PATH=/home/me/vtk_build ..

make

./VectorFieldNonZeroExtraction