<commit_before>/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    FEMImageRegLMEx.cxx
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

  Copyright (c) 2002 Insight Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for detail.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/

// disable debug warnings in MS compiler
#ifdef _MSC_VER
#pragma warning(disable: 4786)
#endif
#include <fstream>
#include "FEMImageRegLMEx.h"

namespace itk {
namespace fem {


/*
 * Generate a rectangular mesh of hexahedron elements
 */
void GenerateRegular3DMesh
(itk::fem::Element::ConstPointer e0, Solver& S, vnl_vector<double>& orig, 
 vnl_vector<double>& size, vnl_vector<double>& Nel)
{

  // Check for correct number of dimensions
  if(orig.size() != Element3DC0LinearHexahedron::NumberOfSpatialDimensions ||
     size.size() != Element3DC0LinearHexahedron::NumberOfSpatialDimensions ||
     Nel.size()  != Element3DC0LinearHexahedron::NumberOfSpatialDimensions)
  {
    throw FEMException(__FILE__, __LINE__, "GenerateMesh<Element2DC0LinearQuadrilateral>::Rectangular");
  }

  // Number of nodes in each dimension
  Nel[0]=floor(Nel[0]);
  Nel[1]=floor(Nel[1]);
  Nel[2]=floor(Nel[2]);
  unsigned int Ni=static_cast<unsigned int>(Nel[0]);
  unsigned int Nj=static_cast<unsigned int>(Nel[1]);
  unsigned int Nk=static_cast<unsigned int>(Nel[2]);

  // Create nodes
  Node::Pointer n;
  int gn=0; // number of node
  for(unsigned int k=0; k<=Nk; k++)
  {
    for(unsigned int j=0; j<=Nj; j++)
    {
      for(unsigned int i=0; i<=Ni; i++)
      {
        n=new Node(orig[0]+i*size[0]/Nel[0],
                   orig[1]+j*size[1]/Nel[1],
                   orig[2]+k*size[2]/Nel[2]);
        n->GN=gn;
        gn++;
        S.node.push_back(FEMP<Node>(n));
      }
    }
  }

  // Create elements  
  gn=0; // global number of the element
  itk::fem::Element3DC0LinearHexahedron::Pointer e;
  for(unsigned int k=0; k<Nk; k++)
  {
    for(unsigned int j=0; j<Nj; j++)
    {
      for(unsigned int i=0; i<Ni; i++)
      {
        e=dynamic_cast<Element3DC0LinearHexahedron*>(e0->Clone());
        e->SetNode(0,S.node.Find( i+  (Ni+1)*(j  +(Nj+1)*k) ));
        e->SetNode(1,S.node.Find( i+1+(Ni+1)*(j  +(Nj+1)*k) ));
        e->SetNode(2,S.node.Find( i+1+(Ni+1)*(j+1+(Nj+1)*k) ));
        e->SetNode(3,S.node.Find( i+  (Ni+1)*(j+1+(Nj+1)*k) ));
        e->SetNode(4,S.node.Find( i+  (Ni+1)*(j  +(Nj+1)*(k+1)) ));
        e->SetNode(5,S.node.Find( i+1+(Ni+1)*(j  +(Nj+1)*(k+1)) ));
        e->SetNode(6,S.node.Find( i+1+(Ni+1)*(j+1+(Nj+1)*(k+1)) ));
        e->SetNode(7,S.node.Find( i+  (Ni+1)*(j+1+(Nj+1)*(k+1)) ));

        e->GN=gn;
        gn++;
        S.el.push_back(FEMP<Element>(e));
      }
    }
  }

}

template<class TReference,class TTarget, class TElement>
ImageRegLMEx<TReference,TTarget,TElement>::~ImageRegLMEx( )
{
 
}

template<class TReference,class TTarget, class TElement>
ImageRegLMEx<TReference,TTarget,TElement>::ImageRegLMEx( )
{
  m_FileCount=0; 
  m_MeshResolution=16*2;  // determines the 'resolution' of the grid
 
  m_MinE=vnl_huge_val(m_MinE);  

  m_DescentDirection=positive;
  m_E=1.; 
  m_Maxiters=1;  
  m_dT=1;
  m_Rho=1;

  m_NumberOfIntegrationPoints=4;
  m_MetricWidth=3;
  m_DoLineSearchOnImageEnergy=false;
  m_SearchForMinAtEachLevel=false;
  m_LineSearchStep=0.025;

  m_NumLevels=1;
  m_MaxLevel=1;
  m_MeshStep=2;
  m_MeshLevels=1;
  m_DoMultiRes=false;
  m_UseLandmarks=false;

  m_ReferenceFileName = NULL;
  m_TargetFileName = NULL;
  m_LandmarkFileName = NULL;
  m_Field=NULL;
  m_TotalIterations=0;

  for (unsigned int i=0; i < ImageDimension; i++)
  {
     m_ImageScaling[i]=1;
     m_ImageSize[i]=0;
     m_ImageOrigin[i]=0;
  }
  
  
  // Register the correct element type and load implementation with the visitor dispatcher 
  VisitorDispatcher<ElementType,itk::fem::Element::LoadType, ElementType::LoadImplementationFunctionPointer>
  ::RegisterVisitor((ImageMetricLoadType*)0, &ImageMetricLoadImplementation<ImageMetricLoadType>::ImplementImageMetricLoad);


}


template<class TReference,class TTarget, class TElement>
void ImageRegLMEx<TReference,TTarget,TElement>::RunRegistration()
{
    std::cout << "beginning \n";  
    m_Solver.SetDeltatT(m_dT);  
    m_Solver.SetRho(m_Rho);    

    ReadImages();
 
   
  // Solve the system in time 

  if (!m_DoMultiRes) 
  {
    
    CreateMesh(m_ImageOrigin,m_ImageSize,(double)m_MeshResolution,m_Solver); 
    ApplyLoads(m_Solver,m_MeshResolution);

    m_Solver.GenerateGFN(); 

    LinearSystemWrapperItpack itpackWrapper; 
    itpackWrapper.SetMaximumNonZeroValuesInMatrix(25*m_Solver.GetNGFN());
    itpackWrapper.SetMaximumNumberIterations(m_Solver.GetNGFN()); 
    itpackWrapper.SetTolerance(1.e-6);
    itpackWrapper.JacobianSemiIterative(); 
    m_Solver.SetLinearSystemWrapper(&itpackWrapper); 

    m_Load=ImageRegLMEx<TReference,TTarget,TElement>::ImageMetricLoadType::New();

    m_Load->SetReferenceImage(m_RefImg);
    m_Load->SetTargetImage(m_TarImg);
    m_Load->SetMetric(m_Metric);
    m_Load->InitializeMetric();
    typename ImageType::SizeType r={{m_MetricWidth,m_MetricWidth}};
    m_Load->SetMetricRadius(r);
    m_Load->SetNumberOfIntegrationPoints(m_NumberOfIntegrationPoints);
    m_Load->GN=m_Solver.load.size()+1; //NOTE SETTING GN FOR FIND LATER
    m_Load->SetSign((Float)m_DescentDirection);
    m_Solver.load.push_back( FEMP<Load>(&*m_Load) );    
    m_Load=dynamic_cast<ImageRegLMEx<TReference,TTarget,TElement>::ImageMetricLoadType*> (&*m_Solver.load.Find(m_Solver.load.size()));  
  
 
    m_Solver.AssembleKandM();
 
    IterativeSolve(m_Solver);

    GetVectorField(m_Solver);
    WarpImage(); 
  }
  else 
  {
    MultiResSolve();
  }
    

  std::cout<<"\n E " << m_E << " dt " << m_dT << " rho " << m_Rho << "\n";
  
}

template<class TReference,class TTarget, class TElement>
void ImageRegLMEx<TReference,TTarget,TElement>::ReadImages()
{
  // DEFINE INPUT FILES
 
  // Read a Raw File
  typedef  itk::ImageFileReader< ImageType >      FileSourceType;
  typedef  itk::RawImageIO< ImageType::PixelType,ImageType::ImageDimension>   RawReaderType;


  FileSourceType::Pointer reffilter = FileSourceType::New();
  reffilter->SetFileName( m_ReferenceFileName );
  FileSourceType::Pointer tarfilter = FileSourceType::New();
  tarfilter->SetFileName( m_TargetFileName );

  RawReaderType::Pointer  rawReader  = RawReaderType::New();
  
  for (int ii=0; ii<ImageDimension; ii++) rawReader->SetDimensions( ii, m_ImageSize[ii] );

  reffilter->SetImageIO( rawReader );
  tarfilter->SetImageIO( rawReader );

  try
    {
    reffilter->Update();
    }
  catch( itk::ExceptionObject & e )
    {
    std::cerr << "Exception caught during reference file reading " << std::endl;
    std::cerr << e << std::endl;
    return;
    }<commit_msg>FIX: missing typenames.<commit_after>/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    FEMImageRegLMEx.cxx
  Language:  C++
  Date:      $Date$
  Version:   $Revision$

  Copyright (c) 2002 Insight Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for detail.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/

// disable debug warnings in MS compiler
#ifdef _MSC_VER
#pragma warning(disable: 4786)
#endif
#include <fstream>
#include "FEMImageRegLMEx.h"

namespace itk {
namespace fem {


/*
 * Generate a rectangular mesh of hexahedron elements
 */
void GenerateRegular3DMesh
(itk::fem::Element::ConstPointer e0, Solver& S, vnl_vector<double>& orig, 
 vnl_vector<double>& size, vnl_vector<double>& Nel)
{

  // Check for correct number of dimensions
  if(orig.size() != Element3DC0LinearHexahedron::NumberOfSpatialDimensions ||
     size.size() != Element3DC0LinearHexahedron::NumberOfSpatialDimensions ||
     Nel.size()  != Element3DC0LinearHexahedron::NumberOfSpatialDimensions)
  {
    throw FEMException(__FILE__, __LINE__, "GenerateMesh<Element2DC0LinearQuadrilateral>::Rectangular");
  }

  // Number of nodes in each dimension
  Nel[0]=floor(Nel[0]);
  Nel[1]=floor(Nel[1]);
  Nel[2]=floor(Nel[2]);
  unsigned int Ni=static_cast<unsigned int>(Nel[0]);
  unsigned int Nj=static_cast<unsigned int>(Nel[1]);
  unsigned int Nk=static_cast<unsigned int>(Nel[2]);

  // Create nodes
  Node::Pointer n;
  int gn=0; // number of node
  for(unsigned int k=0; k<=Nk; k++)
  {
    for(unsigned int j=0; j<=Nj; j++)
    {
      for(unsigned int i=0; i<=Ni; i++)
      {
        n=new Node(orig[0]+i*size[0]/Nel[0],
                   orig[1]+j*size[1]/Nel[1],
                   orig[2]+k*size[2]/Nel[2]);
        n->GN=gn;
        gn++;
        S.node.push_back(FEMP<Node>(n));
      }
    }
  }

  // Create elements  
  gn=0; // global number of the element
  itk::fem::Element3DC0LinearHexahedron::Pointer e;
  for(unsigned int k=0; k<Nk; k++)
  {
    for(unsigned int j=0; j<Nj; j++)
    {
      for(unsigned int i=0; i<Ni; i++)
      {
        e=dynamic_cast<Element3DC0LinearHexahedron*>(e0->Clone());
        e->SetNode(0,S.node.Find( i+  (Ni+1)*(j  +(Nj+1)*k) ));
        e->SetNode(1,S.node.Find( i+1+(Ni+1)*(j  +(Nj+1)*k) ));
        e->SetNode(2,S.node.Find( i+1+(Ni+1)*(j+1+(Nj+1)*k) ));
        e->SetNode(3,S.node.Find( i+  (Ni+1)*(j+1+(Nj+1)*k) ));
        e->SetNode(4,S.node.Find( i+  (Ni+1)*(j  +(Nj+1)*(k+1)) ));
        e->SetNode(5,S.node.Find( i+1+(Ni+1)*(j  +(Nj+1)*(k+1)) ));
        e->SetNode(6,S.node.Find( i+1+(Ni+1)*(j+1+(Nj+1)*(k+1)) ));
        e->SetNode(7,S.node.Find( i+  (Ni+1)*(j+1+(Nj+1)*(k+1)) ));

        e->GN=gn;
        gn++;
        S.el.push_back(FEMP<Element>(e));
      }
    }
  }

}

template<class TReference,class TTarget, class TElement>
ImageRegLMEx<TReference,TTarget,TElement>::~ImageRegLMEx( )
{
 
}

template<class TReference,class TTarget, class TElement>
ImageRegLMEx<TReference,TTarget,TElement>::ImageRegLMEx( )
{
  m_FileCount=0; 
  m_MeshResolution=16*2;  // determines the 'resolution' of the grid
 
  m_MinE=vnl_huge_val(m_MinE);  

  m_DescentDirection=positive;
  m_E=1.; 
  m_Maxiters=1;  
  m_dT=1;
  m_Rho=1;

  m_NumberOfIntegrationPoints=4;
  m_MetricWidth=3;
  m_DoLineSearchOnImageEnergy=false;
  m_SearchForMinAtEachLevel=false;
  m_LineSearchStep=0.025;

  m_NumLevels=1;
  m_MaxLevel=1;
  m_MeshStep=2;
  m_MeshLevels=1;
  m_DoMultiRes=false;
  m_UseLandmarks=false;

  m_ReferenceFileName = NULL;
  m_TargetFileName = NULL;
  m_LandmarkFileName = NULL;
  m_Field=NULL;
  m_TotalIterations=0;

  for (unsigned int i=0; i < ImageDimension; i++)
  {
     m_ImageScaling[i]=1;
     m_ImageSize[i]=0;
     m_ImageOrigin[i]=0;
  }
  
  
  // Register the correct element type and load implementation with the visitor dispatcher 
  typedef VisitorDispatcher< ElementType,
                             itk::fem::Element::LoadType, 
                             typename ElementType::LoadImplementationFunctionPointer >   DispatcherType;
  
  DispatcherType::RegisterVisitor((ImageMetricLoadType*)0, 
                    &ImageMetricLoadImplementation<ImageMetricLoadType>::ImplementImageMetricLoad);


}


template<class TReference,class TTarget, class TElement>
void ImageRegLMEx<TReference,TTarget,TElement>::RunRegistration()
{
    std::cout << "beginning \n";  
    m_Solver.SetDeltatT(m_dT);  
    m_Solver.SetRho(m_Rho);    

    ReadImages();
 
   
  // Solve the system in time 

  if (!m_DoMultiRes) 
  {
    
    CreateMesh(m_ImageOrigin,m_ImageSize,(double)m_MeshResolution,m_Solver); 
    ApplyLoads(m_Solver,m_MeshResolution);

    m_Solver.GenerateGFN(); 

    LinearSystemWrapperItpack itpackWrapper; 
    itpackWrapper.SetMaximumNonZeroValuesInMatrix(25*m_Solver.GetNGFN());
    itpackWrapper.SetMaximumNumberIterations(m_Solver.GetNGFN()); 
    itpackWrapper.SetTolerance(1.e-6);
    itpackWrapper.JacobianSemiIterative(); 
    m_Solver.SetLinearSystemWrapper(&itpackWrapper); 

    m_Load=ImageRegLMEx<TReference,TTarget,TElement>::ImageMetricLoadType::New();

    m_Load->SetReferenceImage(m_RefImg);
    m_Load->SetTargetImage(m_TarImg);
    m_Load->SetMetric(m_Metric);
    m_Load->InitializeMetric();
    typename ImageType::SizeType r={{m_MetricWidth,m_MetricWidth}};
    m_Load->SetMetricRadius(r);
    m_Load->SetNumberOfIntegrationPoints(m_NumberOfIntegrationPoints);
    m_Load->GN=m_Solver.load.size()+1; //NOTE SETTING GN FOR FIND LATER
    m_Load->SetSign((Float)m_DescentDirection);
    m_Solver.load.push_back( FEMP<Load>(&*m_Load) );    
    m_Load=dynamic_cast<ImageRegLMEx<TReference,TTarget,TElement>::ImageMetricLoadType*> (&*m_Solver.load.Find(m_Solver.load.size()));  
  
 
    m_Solver.AssembleKandM();
 
    IterativeSolve(m_Solver);

    GetVectorField(m_Solver);
    WarpImage(); 
  }
  else 
  {
    MultiResSolve();
  }
    

  std::cout<<"\n E " << m_E << " dt " << m_dT << " rho " << m_Rho << "\n";
  
}

template<class TReference,class TTarget, class TElement>
void ImageRegLMEx<TReference,TTarget,TElement>::ReadImages()
{
  // DEFINE INPUT FILES
 
  // Read a Raw File
  typedef  itk::ImageFileReader< ImageType >      FileSourceType;
  typedef  itk::RawImageIO< typename ImageType::PixelType,
                            ImageType::ImageDimension>   RawReaderType;


  FileSourceType::Pointer reffilter = FileSourceType::New();
  reffilter->SetFileName( m_ReferenceFileName );
  FileSourceType::Pointer tarfilter = FileSourceType::New();
  tarfilter->SetFileName( m_TargetFileName );

  RawReaderType::Pointer  rawReader  = RawReaderType::New();
  
  for (int ii=0; ii<ImageDimension; ii++) rawReader->SetDimensions( ii, m_ImageSize[ii] );

  reffilter->SetImageIO( rawReader );
  tarfilter->SetImageIO( rawReader );

  try
    {
    reffilter->Update();
    }
  catch( itk::ExceptionObject & e )
    {
    std::cerr << "Exception caught during reference file reading " << std::endl;
    std::cerr << e << std::endl;
    return;
    }
{'fraction_non_alphanumeric': 0.16249223763196025, 'fraction_numerical': 0.01586973021458635, 'mean_word_length': 4.030892051371052, 'pattern_counts': {'":': 0, '<': 121, '<?xml version=': 0, '>': 107, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '477894', 'n_tokens_mistral': 5433, 'n_tokens_neox': 5183, 'n_words': 843}