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

  Program:   Advanced Normalization Tools
  Module:    $RCSfile: MeasureImageSimilarity.cxx,v $
  Language:  C++
  Date:      $Date: 2009/01/05 20:09:47 $
  Version:   $Revision: 1.19 $

  Copyright (c) ConsortiumOfANTS. All rights reserved.
  See accompanying COPYING.txt or
 http://sourceforge.net/projects/advants/files/ANTS/ANTSCopyright.txt for

     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.

=========================================================================*/
#include "ReadWriteImage.h"
#include "itkDiscreteGaussianImageFilter.h"
#include "itkAvantsMutualInformationRegistrationFunction.h"
#include "itkSpatialMutualInformationRegistrationFunction.h"
#include "itkProbabilisticRegistrationFunction.h"
#include "itkCrossCorrelationRegistrationFunction.h"
// #include "itkLandmarkCrossCorrelationRegistrationFunction.h"

template <unsigned int ImageDimension>
int MeasureImageSimilarity(unsigned int argc, char *argv[])
{
  typedef float                                                  PixelType;
  typedef itk::Vector<float, ImageDimension>                     VectorType;
  typedef itk::Image<VectorType, ImageDimension>                 FieldType;
  typedef itk::Image<PixelType, ImageDimension>                  ImageType;
  typedef itk::ImageFileWriter<ImageType>                        writertype;
  typedef typename ImageType::IndexType                          IndexType;
  typedef typename ImageType::SizeType                           SizeType;
  typedef typename ImageType::SpacingType                        SpacingType;
  typedef itk::AffineTransform<double, ImageDimension>           AffineTransformType;
  typedef itk::LinearInterpolateImageFunction<ImageType, double> InterpolatorType1;
  typedef itk::ImageRegionIteratorWithIndex<ImageType>           Iterator;

  typedef itk::Image<float, 2>                JointHistType;
  typedef itk::ImageFileWriter<JointHistType> jhwritertype;

// get command line params
  unsigned int argct = 2;
  unsigned int whichmetric = atoi(argv[argct]); argct++;
  std::string  fn1 = std::string(argv[argct]); argct++;
  std::string  fn2 = std::string(argv[argct]); argct++;
  std::string  logfilename = "";
  if( argc > argct )
    {
    logfilename = std::string(argv[argct]);
    }
  argct++;
  std::string imgfilename = "";
  if( argc > argct )
    {
    imgfilename = std::string(argv[argct]);
    }
  argct++;
  double targetvalue = 0;
  if( argc > argct )
    {
    targetvalue = atof(argv[argct]);
    }
  argct++;
  double epsilontolerance = 1.e20;
  if( argc > argct )
    {
    epsilontolerance = atof(argv[argct]);
    }
  argct++;

  typename ImageType::Pointer image1 = NULL;
  ReadImage<ImageType>(image1, fn1.c_str() );
  typename ImageType::Pointer image2 = NULL;
  ReadImage<ImageType>(image2, fn2.c_str() );

/*
  typedef itk::ImageRegionIteratorWithIndex<FieldType> VIterator;
  typename FieldType::Pointer field=FieldType::New();
  field->SetLargestPossibleRegion( image1->GetLargestPossibleRegion() );
  field->SetBufferedRegion( image1->GetLargestPossibleRegion() );
  field->SetLargestPossibleRegion( image1->GetLargestPossibleRegion() );
  field->Allocate();
  field->SetSpacing(image1->GetSpacing());
  field->SetOrigin(image1->GetOrigin());
  VectorType zero;
  zero.Fill(0);
  VIterator vfIter2( field,  field->GetLargestPossibleRegion() );
  for(  vfIter2.GoToBegin(); !vfIter2.IsAtEnd(); ++vfIter2 )
    {
      IndexType index=vfIter2.GetIndex();
      vfIter2.Set(zero);
    }
*/
  typename ImageType::Pointer metricimg = ImageType::New();
  metricimg->SetRegions( image1->GetLargestPossibleRegion() );
  metricimg->Allocate();
  metricimg->SetSpacing(image1->GetSpacing() );
  metricimg->SetOrigin(image1->GetOrigin() );
  metricimg->SetDirection(image1->GetDirection() );
  Iterator iter( metricimg,  metricimg->GetLargestPossibleRegion() );
  for(  iter.GoToBegin(); !iter.IsAtEnd(); ++iter )
    {
    iter.Set(0);
    }

  typedef ImageType FixedImageType;
  typedef ImageType MovingImageType;
  typedef FieldType DisplacementFieldType;

  // Choose the similarity metric
  typedef itk::AvantsMutualInformationRegistrationFunction<FixedImageType, MovingImageType,
                                                           DisplacementFieldType>  MIMetricType;
  typedef itk::SpatialMutualInformationRegistrationFunction<FixedImageType, MovingImageType,
                                                            DisplacementFieldType> SMIMetricType;
  typedef itk::CrossCorrelationRegistrationFunction<FixedImageType, MovingImageType,
                                                    DisplacementFieldType>         CCMetricType;
  // typedef itk::LandmarkCrossCorrelationRegistrationFunction<FixedImageType,MovingImageType,DisplacementFieldType>
  // MetricType;
  // typename
  typename MIMetricType::Pointer mimet = MIMetricType::New();
  typename SMIMetricType::Pointer smimet = SMIMetricType::New();
  typename CCMetricType::Pointer ccmet = CCMetricType::New();

//  int nbins=32;

  typename CCMetricType::RadiusType hradius;
  typename CCMetricType::RadiusType ccradius;
  ccradius.Fill(4);
  typename MIMetricType::RadiusType miradius;
  miradius.Fill(0);

//  mimet->SetDisplacementField(field);
  mimet->SetFixedImage(image1);
  mimet->SetMovingImage(image2);
  mimet->SetRadius(miradius);
  mimet->SetGradientStep(1.e2);
  mimet->SetNormalizeGradient(false);

//  ccmet->SetDisplacementField(field);
  ccmet->SetFixedImage(image1);
  ccmet->SetMovingImage(image2);
  ccmet->SetRadius(ccradius);
  ccmet->SetGradientStep(1.e2);
  ccmet->SetNormalizeGradient(false);

  double      metricvalue = 0;
  std::string metricname = "";
  if( whichmetric  == 0 )
    {
    hradius = miradius;
    unsigned long ct = 0;
    for(  iter.GoToBegin(); !iter.IsAtEnd(); ++iter )
      {
      IndexType index = iter.GetIndex();
      double    fval = image1->GetPixel(index);
      double    mval = image2->GetPixel(index);
      metricvalue += fabs(fval - mval);
      metricimg->SetPixel(index, fabs(fval - mval) );
      ct++;
      }
    metricvalue /= (float)ct;
    metricname = "MSQ ";
    }
  else if( whichmetric == 1 ) // imagedifference
    {
    hradius = ccradius;
    metricname = "CC ";
    std::cout << metricname << std::endl;
    ccmet->InitializeIteration();
    metricimg = ccmet->MakeImage();
    metricvalue = ccmet->ComputeCrossCorrelation() * (1.0);
    }
  else if( whichmetric == 2 )
    {
    hradius = miradius;
    mimet->InitializeIteration();
    metricvalue = mimet->ComputeMutualInformation();
    metricname = "MI ";
    }
  else if( whichmetric == 3 )
    {
    hradius = miradius;
    smimet->InitializeIteration();
    metricvalue = smimet->ComputeSpatialMutualInformation();
    metricname = "SMI ";
    }
  std::cout << fn1 << " : " << fn2 << " => " <<  metricname << metricvalue << std::endl;
  if( logfilename.length() > 3 )
    {
    std::ofstream logfile;
    logfile.open(logfilename.c_str(), std::ofstream::app);
    if( logfile.good() )
      {
      logfile <<  fn1 << " : " << fn2 << " => " << metricname << metricvalue << std::endl;<commit_msg>ENH: cc image metric update <commit_after>/*=========================================================================

  Program:   Advanced Normalization Tools
  Module:    $RCSfile: MeasureImageSimilarity.cxx,v $
  Language:  C++
  Date:      $Date: 2009/01/05 20:09:47 $
  Version:   $Revision: 1.19 $

  Copyright (c) ConsortiumOfANTS. All rights reserved.
  See accompanying COPYING.txt or
 http://sourceforge.net/projects/advants/files/ANTS/ANTSCopyright.txt for

     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.

=========================================================================*/
#include "ReadWriteImage.h"
#include "itkDiscreteGaussianImageFilter.h"
#include "itkAvantsMutualInformationRegistrationFunction.h"
#include "itkSpatialMutualInformationRegistrationFunction.h"
#include "itkProbabilisticRegistrationFunction.h"
#include "itkCrossCorrelationRegistrationFunction.h"
// #include "itkLandmarkCrossCorrelationRegistrationFunction.h"

template <unsigned int ImageDimension>
int MeasureImageSimilarity(unsigned int argc, char *argv[])
{
  typedef float                                                  PixelType;
  typedef itk::Vector<float, ImageDimension>                     VectorType;
  typedef itk::Image<VectorType, ImageDimension>                 FieldType;
  typedef itk::Image<PixelType, ImageDimension>                  ImageType;
  typedef itk::ImageFileWriter<ImageType>                        writertype;
  typedef typename ImageType::IndexType                          IndexType;
  typedef typename ImageType::SizeType                           SizeType;
  typedef typename ImageType::SpacingType                        SpacingType;
  typedef itk::AffineTransform<double, ImageDimension>           AffineTransformType;
  typedef itk::LinearInterpolateImageFunction<ImageType, double> InterpolatorType1;
  typedef itk::ImageRegionIteratorWithIndex<ImageType>           Iterator;

  typedef itk::Image<float, 2>                JointHistType;
  typedef itk::ImageFileWriter<JointHistType> jhwritertype;

// get command line params
  unsigned int argct = 2;
  unsigned int whichmetric = atoi(argv[argct]); argct++;
  std::string  fn1 = std::string(argv[argct]); argct++;
  std::string  fn2 = std::string(argv[argct]); argct++;
  std::string  logfilename = "";
  if( argc > argct )
    {
    logfilename = std::string(argv[argct]);
    }
  argct++;
  std::string imgfilename = "";
  if( argc > argct )
    {
    imgfilename = std::string(argv[argct]);
    }
  argct++;
  double targetvalue = 0;
  if( argc > argct )
    {
    targetvalue = atof(argv[argct]);
    }
  argct++;
  double epsilontolerance = 1.e20;
  if( argc > argct )
    {
    epsilontolerance = atof(argv[argct]);
    }
  argct++;

  typename ImageType::Pointer image1 = NULL;
  ReadImage<ImageType>(image1, fn1.c_str() );
  typename ImageType::Pointer image2 = NULL;
  ReadImage<ImageType>(image2, fn2.c_str() );

/*
  typedef itk::ImageRegionIteratorWithIndex<FieldType> VIterator;
  typename FieldType::Pointer field=FieldType::New();
  field->SetLargestPossibleRegion( image1->GetLargestPossibleRegion() );
  field->SetBufferedRegion( image1->GetLargestPossibleRegion() );
  field->SetLargestPossibleRegion( image1->GetLargestPossibleRegion() );
  field->Allocate();
  field->SetSpacing(image1->GetSpacing());
  field->SetOrigin(image1->GetOrigin());
  VectorType zero;
  zero.Fill(0);
  VIterator vfIter2( field,  field->GetLargestPossibleRegion() );
  for(  vfIter2.GoToBegin(); !vfIter2.IsAtEnd(); ++vfIter2 )
    {
      IndexType index=vfIter2.GetIndex();
      vfIter2.Set(zero);
    }
*/
  typename ImageType::Pointer metricimg = ImageType::New();
  metricimg->SetRegions( image1->GetLargestPossibleRegion() );
  metricimg->Allocate();
  metricimg->SetSpacing(image1->GetSpacing() );
  metricimg->SetOrigin(image1->GetOrigin() );
  metricimg->SetDirection(image1->GetDirection() );
  Iterator iter( metricimg,  metricimg->GetLargestPossibleRegion() );
  for(  iter.GoToBegin(); !iter.IsAtEnd(); ++iter )
    {
    iter.Set(0);
    }

  typedef ImageType FixedImageType;
  typedef ImageType MovingImageType;
  typedef FieldType DisplacementFieldType;

  // Choose the similarity metric
  typedef itk::AvantsMutualInformationRegistrationFunction<FixedImageType, MovingImageType,
                                                           DisplacementFieldType>  MIMetricType;
  typedef itk::SpatialMutualInformationRegistrationFunction<FixedImageType, MovingImageType,
                                                            DisplacementFieldType> SMIMetricType;
  typedef itk::CrossCorrelationRegistrationFunction<FixedImageType, MovingImageType,
                                                    DisplacementFieldType>         CCMetricType;
  // typedef itk::LandmarkCrossCorrelationRegistrationFunction<FixedImageType,MovingImageType,DisplacementFieldType>
  // MetricType;
  // typename
  typename MIMetricType::Pointer mimet = MIMetricType::New();
  typename SMIMetricType::Pointer smimet = SMIMetricType::New();
  typename CCMetricType::Pointer ccmet = CCMetricType::New();

//  int nbins=32;

  typename CCMetricType::RadiusType hradius;
  typename CCMetricType::RadiusType ccradius;
  ccradius.Fill(4);
  typename MIMetricType::RadiusType miradius;
  miradius.Fill(0);

//  mimet->SetDisplacementField(field);
  mimet->SetFixedImage(image1);
  mimet->SetMovingImage(image2);
  mimet->SetRadius(miradius);
  mimet->SetGradientStep(1.e2);
  mimet->SetNormalizeGradient(false);

//  ccmet->SetDisplacementField(field);
  ccmet->SetFixedImage(image1);
  ccmet->SetMovingImage(image2);
  ccmet->SetRadius(ccradius);
  ccmet->SetGradientStep(1.e2);
  ccmet->SetNormalizeGradient(false);

  double      metricvalue = 0;
  std::string metricname = "";
  if( whichmetric  == 0 )
    {
    hradius = miradius;
    unsigned long ct = 0;
    for(  iter.GoToBegin(); !iter.IsAtEnd(); ++iter )
      {
      IndexType index = iter.GetIndex();
      double    fval = image1->GetPixel(index);
      double    mval = image2->GetPixel(index);
      metricvalue += fabs(fval - mval);
      metricimg->SetPixel(index, fabs(fval - mval) );
      ct++;
      }
    metricvalue /= (float)ct;
    metricname = "MSQ ";
    }
  else if( whichmetric == 1 ) // imagedifference
    {
    double ccval = 0;
    hradius = ccradius;
    metricname = "CC ";
    typedef itk::ConstNeighborhoodIterator<FixedImageType> ScanIteratorType;
    typename FixedImageType::RegionType region = image1->GetLargestPossibleRegion();
    ScanIteratorType   asamIt( hradius, image1, region);
    unsigned long      ct = 0;
    const unsigned int hoodlen = asamIt.Size();
    double             totval = 0;
    for(  iter.GoToBegin(); !iter.IsAtEnd(); ++iter )
      {
      IndexType index = iter.GetIndex();
      double    val = 0;
      double    fmip = 0, mmip = 0, ffip = 0;
      asamIt.SetLocation(index);
      for( unsigned int i = 0; i < asamIt.Size(); i++ )
        {
        IndexType locind = asamIt.GetIndex(i);
        if( region.IsInside( locind ) )
          {
          double f = image1->GetPixel(locind);
          double m = image2->GetPixel(locind);
          fmip += (f * m);  mmip += (m * m);  ffip += (f * f);
          }
        }
      double denom = mmip * ffip;
      if( denom == 0 )
        {
        val = 1;
        }
      else
        {
        val = fmip / sqrt(denom);
        }
      metricimg->SetPixel(index, val);
      // if (ct % 10000 == 0)
      //        std::cout << val << " index " << index << std::endl;
      //      asamIt.SetLocation(index);
      //      totval+=met->localProbabilistic;
      ccval += val;
      ct++;
      }
    if( ct >  0 )
      {
      metricvalue = ccval / ct;
      }
    else
      {
      metricvalue = 0;
      }
    /*
    std::cout << metricname << std::endl;
    ccmet->InitializeIteration();
    metricimg=ccmet->MakeImage();
    metricvalue=ccmet->ComputeCrossCorrelation()*(1.0);
    */
    }
  else if( whichmetric == 2 )
    {
    hradius = miradius;
    mimet->InitializeIteration();
    metricvalue = mimet->ComputeMutualInformation();
    metricname = "MI ";
    }
  else if( whichmetric == 3 )
    {
    hradius = miradius;
    smimet->InitializeIteration();
    metricvalue = smimet->ComputeSpatialMutualInformation();
    metricname = "SMI ";
    }
  std::cout << fn1 << " : " << fn2 << " => " <<  metricname << metricvalue << std::endl;
  if( logfilename.length() > 3 )
    {
    std::ofstream logfile;
    logfile.open(logfilename.c_str(), std::ofstream::app);
    if( logfile.good() )
      {
      logfile <<  fn1 << " : " << fn2 << " => " << metricname << metricvalue << std::endl;
{'fraction_non_alphanumeric': 0.13744896696771, 'fraction_numerical': 0.010763330446616355, 'mean_word_length': 3.2433070866141733, 'pattern_counts': {'":': 0, '<': 111, '<?xml version=': 0, '>': 146, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '651396', 'n_tokens_mistral': 5053, 'n_tokens_neox': 4827, 'n_words': 937}