<commit_before>#include "antsCommandLineParser.h"
#include "antsUtilities.h"
#include "ReadWriteImage.h"

#include "itkAffineTransform.h"
#include "itkAntiAliasBinaryImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkImageToVTKImageFilter.h"

#include "vtkExtractEdges.h"
#include "vtkMarchingCubes.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkPolyDataConnectivityFilter.h"
#include "vtkPolyDataNormals.h"
#include "vtkSmartPointer.h"
#include "vtkUnsignedCharArray.h"
#include "vtkWindowedSincPolyDataFilter.h"
#include "vtkPolyDataWriter.h"
#include "vtkSTLWriter.h"

#include <vector>

namespace ants
{

float CalculateGenus( vtkPolyData *mesh, bool verbose )
{
  vtkSmartPointer<vtkPolyDataConnectivityFilter> connectivityFilter =
    vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
  connectivityFilter->SetExtractionModeToLargestRegion();
  connectivityFilter->SetInputData( mesh );
  connectivityFilter->Update();

  vtkSmartPointer<vtkExtractEdges> edgeExtractor = vtkSmartPointer<vtkExtractEdges>::New();
  edgeExtractor->SetInputData( connectivityFilter->GetOutput() );
  edgeExtractor->Update();
  vtkPolyData* edges = edgeExtractor->GetOutput();

  vtkIdType numberOfEdges = mesh->GetNumberOfCells();
  vtkIdType numberOfVertices = mesh->GetNumberOfPoints();
  vtkIdType numberOfFaces = mesh->GetNumberOfPolys();

  float genus = 0.5 * ( 2.0 - numberOfVertices + numberOfEdges - numberOfFaces );

  if( verbose )
    {
    std::cout << "Genus = "  << genus << std::endl;
    std::cout << "  number of vertices = "  << numberOfVertices << std::endl;
    std::cout << "  number of edges = "  << numberOfEdges << std::endl;
    std::cout << "  number of faces = "  << numberOfFaces << std::endl;
    }

  return genus;
}

int antsSurfaceFunction( itk::ants::CommandLineParser *parser )
{
  const unsigned int ImageDimension = 3;

  typedef float                                      RealType;
  typedef itk::Image<RealType, ImageDimension>       ImageType;
  typedef itk::Image<int, ImageDimension>            MaskImageType;

  typedef unsigned char                              RgbComponentType;
  typedef itk::RGBPixel<RgbComponentType>            RgbPixelType;
  typedef itk::Image<RgbPixelType, ImageDimension>   RgbImageType;

  // Read in input surface image

  ImageType::Pointer inputImage = NULL;

  RealType defaultColorRed = 255.0;
  RealType defaultColorGreen = 255.0;
  RealType defaultColorBlue = 255.0;

  itk::ants::CommandLineParser::OptionType::Pointer inputImageOption =
    parser->GetOption( "surface-image" );
  if( inputImageOption && inputImageOption->GetNumberOfFunctions() )
    {
    if( inputImageOption->GetFunction( 0 )->GetNumberOfParameters() == 0 )
      {
      std::string inputFile = inputImageOption->GetFunction( 0 )->GetName();
      ReadImage<ImageType>( inputImage, inputFile.c_str() );
      }
    else
      {
      std::string inputFile = inputImageOption->GetFunction( 0 )->GetParameter( 0 );
      ReadImage<ImageType>( inputImage, inputFile.c_str() );

      if( inputImageOption->GetFunction( 0 )->GetNumberOfParameters() > 1 )
        {
        std::vector<RealType> defaultColors = parser->ConvertVector<RealType>(
          inputImageOption->GetFunction( 0 )->GetParameter( 1 ) );
        if( defaultColors.size() == 1 )
          {
          defaultColorRed = defaultColors[0];
          defaultColorGreen = defaultColors[0];
          defaultColorBlue = defaultColors[0];
          }
        else if( defaultColors.size() == 3 )
          {
          defaultColorRed = defaultColors[0];
          defaultColorGreen = defaultColors[1];
          defaultColorBlue = defaultColors[2];
          }
        else
          {
          std::cerr << "Incorrect color format specified." << std::endl;
          return EXIT_FAILURE;
          }
        }
      }
    }
  else
    {
    std::cerr << "Input image not specified." << std::endl;
    return EXIT_FAILURE;
    }


  // There's a reorientation issue between itk image physical space and the mesh space
  // for which we have to account.  See
  // http://www.vtk.org/pipermail/vtkusers/2011-July/068595.html
  //   and
  // http://www.vtk.org/Wiki/VTK/ExamplesBoneYard/Cxx/VolumeRendering/itkVtkImageConvert

  typedef itk::AffineTransform<RealType> RigidTransformType;
  RigidTransformType::Pointer meshToItkImageTransform = RigidTransformType::New();
  RigidTransformType::OutputVectorType offset;
  offset[0] = -inputImage->GetOrigin()[0];
  offset[1] = -inputImage->GetOrigin()[1];
  offset[2] = -inputImage->GetOrigin()[2];
  RigidTransformType::MatrixType matrix;
  for( unsigned int i = 0; i < ImageDimension; i++ )
    {
    for( unsigned int j = 0; j < ImageDimension; j++ )
      {
      matrix( i, j ) = inputImage->GetDirection()( i, j );
      }
    }
  meshToItkImageTransform->SetMatrix( matrix );
  meshToItkImageTransform->SetOffset( offset );

  // Get anti-alias RMSE parameter

  RealType antiAliasRmseParameter = 0.03;

  itk::ants::CommandLineParser::OptionType::Pointer antiAliasRmseOption =
    parser->GetOption( "anti-alias-rmse" );
  if( antiAliasRmseOption && antiAliasRmseOption->GetNumberOfFunctions() )
    {
    antiAliasRmseParameter = parser->Convert<RealType>( antiAliasRmseOption->GetFunction( 0 )->GetName() );
    }

  typedef itk::AntiAliasBinaryImageFilter<ImageType, ImageType> AntiAliasFilterType;
  AntiAliasFilterType::Pointer antiAlias = AntiAliasFilterType::New();
  antiAlias->SetMaximumRMSError( antiAliasRmseParameter );
  antiAlias->SetInput( inputImage );
  antiAlias->Update();

  // Reconstruct binary surface.

  typedef itk::ImageToVTKImageFilter<ImageType> ConnectorType;
  ConnectorType::Pointer connector = ConnectorType::New();
  connector->SetInput( antiAlias->GetOutput() );
  connector->Update();

  vtkSmartPointer<vtkMarchingCubes> marchingCubes = vtkSmartPointer<vtkMarchingCubes>::New();
  marchingCubes->SetInputData( connector->GetOutput() );
  marchingCubes->ComputeScalarsOff();
  marchingCubes->ComputeGradientsOff();
  marchingCubes->SetNumberOfContours( 1 );
  marchingCubes->SetValue( 0, 0.0 );
  marchingCubes->Update();

  // Smooth mesh

  vtkSmartPointer<vtkWindowedSincPolyDataFilter> meshSmoother =
    vtkSmartPointer<vtkWindowedSincPolyDataFilter>::New();
  meshSmoother->SetInputData( marchingCubes->GetOutput() );
  meshSmoother->SetNumberOfIterations( 25 );
  meshSmoother->BoundarySmoothingOff();
  meshSmoother->FeatureEdgeSmoothingOff();
  meshSmoother->SetFeatureAngle( 120.0 );
  meshSmoother->SetPassBand( 0.01 );
  meshSmoother->NonManifoldSmoothingOn();
  meshSmoother->NormalizeCoordinatesOn();
  meshSmoother->Update();

  vtkPolyData *vtkMesh = meshSmoother->GetOutput();

  CalculateGenus( vtkMesh, true );

  // Add the functional overlays

  std::vector<RgbImageType::Pointer>    functionalRgbImages;
  std::vector<MaskImageType::Pointer>   functionalMaskImages;
  std::vector<RealType>                 functionalAlphaValues;

  itk::ants::CommandLineParser::OptionType::Pointer functionalOverlayOption =
    parser->GetOption( "functional-overlay" );
  if( functionalOverlayOption && functionalOverlayOption->GetNumberOfFunctions() )
    {
    for( unsigned int n = 0; n < functionalOverlayOption->GetNumberOfFunctions(); n++ )
      {
      if( functionalOverlayOption->GetFunction( n )->GetNumberOfParameters() < 2 )
        {
        std::cerr << "Error:  each functional overlay must have an RGB image and mask."<commit_msg>WIP:  Getting closer to real genus.<commit_after>#include "antsCommandLineParser.h"
#include "antsUtilities.h"
#include "ReadWriteImage.h"

#include "itkAffineTransform.h"
#include "itkAntiAliasBinaryImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkImageToVTKImageFilter.h"

#include "vtkExtractEdges.h"
#include "vtkMarchingCubes.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkPolyDataConnectivityFilter.h"
#include "vtkPolyDataNormals.h"
#include "vtkSmartPointer.h"
#include "vtkTriangleFilter.h"
#include "vtkUnsignedCharArray.h"
#include "vtkWindowedSincPolyDataFilter.h"
#include "vtkPolyDataWriter.h"
#include "vtkSTLWriter.h"

#include <vector>

namespace ants
{

float CalculateGenus( vtkPolyData *mesh, bool verbose )
{
  vtkSmartPointer<vtkExtractEdges> extractEdges = vtkSmartPointer<vtkExtractEdges>::New();
  extractEdges->SetInputData( mesh );
  extractEdges->Update();

  float numberOfEdges = static_cast<float>( extractEdges->GetOutput()->GetNumberOfLines() );
  float numberOfVertices = static_cast<float>( mesh->GetNumberOfPoints() );
  float numberOfFaces = static_cast<float>( mesh->GetNumberOfPolys() );

  float genus = 0.5 * ( 2.0 - numberOfVertices + numberOfEdges - numberOfFaces );

  if( verbose )
    {
    std::cout << "Genus = "  << genus << std::endl;
    std::cout << "  number of vertices = "  << numberOfVertices << std::endl;
    std::cout << "  number of edges = "  << numberOfEdges << std::endl;
    std::cout << "  number of faces = "  << numberOfFaces << std::endl;
    }

  return genus;
}

int antsSurfaceFunction( itk::ants::CommandLineParser *parser )
{
  const unsigned int ImageDimension = 3;

  typedef float                                      RealType;
  typedef itk::Image<RealType, ImageDimension>       ImageType;
  typedef itk::Image<int, ImageDimension>            MaskImageType;

  typedef unsigned char                              RgbComponentType;
  typedef itk::RGBPixel<RgbComponentType>            RgbPixelType;
  typedef itk::Image<RgbPixelType, ImageDimension>   RgbImageType;

  // Read in input surface image

  ImageType::Pointer inputImage = NULL;

  RealType defaultColorRed = 255.0;
  RealType defaultColorGreen = 255.0;
  RealType defaultColorBlue = 255.0;

  itk::ants::CommandLineParser::OptionType::Pointer inputImageOption =
    parser->GetOption( "surface-image" );
  if( inputImageOption && inputImageOption->GetNumberOfFunctions() )
    {
    if( inputImageOption->GetFunction( 0 )->GetNumberOfParameters() == 0 )
      {
      std::string inputFile = inputImageOption->GetFunction( 0 )->GetName();
      ReadImage<ImageType>( inputImage, inputFile.c_str() );
      }
    else
      {
      std::string inputFile = inputImageOption->GetFunction( 0 )->GetParameter( 0 );
      ReadImage<ImageType>( inputImage, inputFile.c_str() );

      if( inputImageOption->GetFunction( 0 )->GetNumberOfParameters() > 1 )
        {
        std::vector<RealType> defaultColors = parser->ConvertVector<RealType>(
          inputImageOption->GetFunction( 0 )->GetParameter( 1 ) );
        if( defaultColors.size() == 1 )
          {
          defaultColorRed = defaultColors[0];
          defaultColorGreen = defaultColors[0];
          defaultColorBlue = defaultColors[0];
          }
        else if( defaultColors.size() == 3 )
          {
          defaultColorRed = defaultColors[0];
          defaultColorGreen = defaultColors[1];
          defaultColorBlue = defaultColors[2];
          }
        else
          {
          std::cerr << "Incorrect color format specified." << std::endl;
          return EXIT_FAILURE;
          }
        }
      }
    }
  else
    {
    std::cerr << "Input image not specified." << std::endl;
    return EXIT_FAILURE;
    }


  // There's a reorientation issue between itk image physical space and the mesh space
  // for which we have to account.  See
  // http://www.vtk.org/pipermail/vtkusers/2011-July/068595.html
  //   and
  // http://www.vtk.org/Wiki/VTK/ExamplesBoneYard/Cxx/VolumeRendering/itkVtkImageConvert

  typedef itk::AffineTransform<RealType> RigidTransformType;
  RigidTransformType::Pointer meshToItkImageTransform = RigidTransformType::New();
  RigidTransformType::OutputVectorType offset;
  offset[0] = -inputImage->GetOrigin()[0];
  offset[1] = -inputImage->GetOrigin()[1];
  offset[2] = -inputImage->GetOrigin()[2];
  RigidTransformType::MatrixType matrix;
  for( unsigned int i = 0; i < ImageDimension; i++ )
    {
    for( unsigned int j = 0; j < ImageDimension; j++ )
      {
      matrix( i, j ) = inputImage->GetDirection()( i, j );
      }
    }
  meshToItkImageTransform->SetMatrix( matrix );
  meshToItkImageTransform->SetOffset( offset );

  // Get anti-alias RMSE parameter

  RealType antiAliasRmseParameter = 0.03;

  itk::ants::CommandLineParser::OptionType::Pointer antiAliasRmseOption =
    parser->GetOption( "anti-alias-rmse" );
  if( antiAliasRmseOption && antiAliasRmseOption->GetNumberOfFunctions() )
    {
    antiAliasRmseParameter = parser->Convert<RealType>( antiAliasRmseOption->GetFunction( 0 )->GetName() );
    }

  typedef itk::AntiAliasBinaryImageFilter<ImageType, ImageType> AntiAliasFilterType;
  AntiAliasFilterType::Pointer antiAlias = AntiAliasFilterType::New();
  antiAlias->SetMaximumRMSError( antiAliasRmseParameter );
  antiAlias->SetInput( inputImage );
  antiAlias->Update();

  // Reconstruct binary surface.

  typedef itk::ImageToVTKImageFilter<ImageType> ConnectorType;
  ConnectorType::Pointer connector = ConnectorType::New();
  connector->SetInput( antiAlias->GetOutput() );
  connector->Update();

  vtkSmartPointer<vtkMarchingCubes> marchingCubes = vtkSmartPointer<vtkMarchingCubes>::New();
  marchingCubes->SetInputData( connector->GetOutput() );
  marchingCubes->ComputeScalarsOff();
  marchingCubes->ComputeGradientsOff();
  marchingCubes->SetNumberOfContours( 1 );
  marchingCubes->SetValue( 0, 0.0 );
  marchingCubes->Update();

  vtkSmartPointer<vtkPolyDataConnectivityFilter> connectivityFilter =
    vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
  connectivityFilter->SetExtractionModeToLargestRegion();
  connectivityFilter->SetInputData( marchingCubes->GetOutput() );
  connectivityFilter->Update();

  vtkSmartPointer<vtkTriangleFilter> triangularizer = vtkSmartPointer<vtkTriangleFilter>::New();
  triangularizer->SetInputData( connectivityFilter->GetOutput() );
  triangularizer->Update();

  // Smooth mesh

  vtkSmartPointer<vtkWindowedSincPolyDataFilter> meshSmoother =
    vtkSmartPointer<vtkWindowedSincPolyDataFilter>::New();
  meshSmoother->SetInputData( triangularizer->GetOutput() );
  meshSmoother->SetNumberOfIterations( 25 );
  meshSmoother->BoundarySmoothingOff();
  meshSmoother->FeatureEdgeSmoothingOff();
  meshSmoother->SetFeatureAngle( 120.0 );
  meshSmoother->SetPassBand( 0.01 );
  meshSmoother->NonManifoldSmoothingOn();
  meshSmoother->NormalizeCoordinatesOn();
  meshSmoother->Update();

  vtkPolyData *vtkMesh = meshSmoother->GetOutput();

  CalculateGenus( vtkMesh, true );

  // Add the functional overlays

  std::vector<RgbImageType::Pointer>    functionalRgbImages;
  std::vector<MaskImageType::Pointer>   functionalMaskImages;
  std::vector<RealType>                 functionalAlphaValues;

  itk::ants::CommandLineParser::OptionType::Pointer functionalOverlayOption =
    parser->GetOption( "functional-overlay" );
  if( functionalOverlayOption && functionalOverlayOption->GetNumberOfFunctions() )
    {
    for( unsigned int n = 0; n < functionalOverlayOption->GetNumberOfFunctions(); n++ )
      {
      if( functionalOverlayOption->GetFunction( n )->GetNumberOfParameters() < 2 )
        {
        std::cerr << "Error:  each functional overlay must have an RGB image and mask."
{'fraction_non_alphanumeric': 0.12178522348013873, 'fraction_numerical': 0.009292585563772004, 'mean_word_length': 4.143722652305621, 'pattern_counts': {'":': 0, '<': 132, '<?xml version=': 0, '>': 187, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4502538', 'n_tokens_mistral': 4922, 'n_tokens_neox': 4713, 'n_words': 852}