<|start_filename|>Tools/DICOMTools/DICOM2NRRDConverterHCC.cxx<|end_filename|>
#include "itkImage.h"
#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkImageSeriesReader.h"
#include "itkImageSeriesWriter.h"

#include "itkAddImageFilter.h"
#include "itkImageFileWriter.h"

#include "itksys/SystemTools.hxx"

#include <boost/algorithm/string.hpp>



int refineString(std::string &in)
{
    boost::trim(in);
    boost::replace_all(in, ",", "");
    boost::replace_all(in, " ", "_");

    return in.length();
}

std::string refineName(std::string &fullName)
{
    std::vector<std::string> strs;
    std::string refineName;

    boost::algorithm::to_lower(fullName);
    boost::split(strs, fullName, boost::is_any_of("^"));

    boost::trim_if(strs[0], ! boost::algorithm::is_alpha());
    boost::trim_if(strs[1], ! boost::algorithm::is_alpha());

    strs[0][0] = toupper(strs[0][0]);
    strs[1][0] = toupper(strs[1][0]);

    refineName = strs[1] + "_" + strs[0];

    return refineName;
}

std::string initialName(std::string &fullName)
{
    std::vector<std::string> strs;
    std::string initialName;

    boost::algorithm::to_lower(fullName);
    boost::split(strs, fullName, boost::is_any_of("^"));

    boost::trim_if(strs[0], ! boost::algorithm::is_alpha());
    boost::trim_if(strs[1], ! boost::algorithm::is_alpha());

    strs[0][0] = toupper(strs[0][0]);
    strs[1][0] = toupper(strs[1][0]);

    initialName = strs[1][0] + strs[0];

    return initialName;
}


int main( int argc, char *argv[] )
{
    if ( argc < 3 )
    {
        std::cerr << "Usage: " << std::endl;
        std::cerr << argv[0] << " DicomDirectory OutputDicomDirectory"
                  << std::endl;
        return EXIT_FAILURE;
    }

    using PixelType = signed short;
    using ImageType = itk::Image< PixelType, 3 >;
    using DICOMImageType = itk::Image< PixelType, 2 >;

    using NamesGeneratorType = itk::GDCMSeriesFileNames;
    NamesGeneratorType::Pointer nameGenerator = NamesGeneratorType::New();
    nameGenerator->SetUseSeriesDetails( true );
    nameGenerator->AddSeriesRestriction("0008|0021" );
    nameGenerator->SetRecursive( true );
    nameGenerator->SetDirectory( argv[1] );

    try
    {
        std::cout << std::endl << "The directory: " << std::endl;
        std::cout << std::endl << argv[1] << std::endl << std::endl;
        std::cout << "Contains the following DICOM Series: ";
        std::cout << std::endl << std::endl;


        using SeriesIdContainer = std::vector< std::string >;
        const SeriesIdContainer &seriesUID = nameGenerator->GetSeriesUIDs();
        SeriesIdContainer::const_iterator seriesItr = seriesUID.begin();
        SeriesIdContainer::const_iterator seriesEnd = seriesUID.end();

        while ( seriesItr != seriesEnd )
        {
            std::cout << seriesItr->c_str() << std::endl;
            ++seriesItr;
        }

        std::string studyDescriptionOld;

        for ( seriesItr = seriesUID.begin(); seriesItr != seriesEnd; seriesItr++)
        {
            std::string seriesIdentifier;
            seriesIdentifier = seriesItr->c_str();

            std::cout << std::endl << std::endl;
            std::cout << "Now reading series: " << std::endl << std::endl;
            std::cout << seriesIdentifier << std::endl;
            std::cout << std::endl << std::endl;

            using ReaderType = itk::ImageSeriesReader< ImageType >;
            ReaderType::Pointer reader = ReaderType::New();

            using ImageIOType = itk::GDCMImageIO;
            ImageIOType::Pointer dicomIO = ImageIOType::New();
            reader->SetImageIO( dicomIO );

            using FileNamesContainer = std::vector< std::string >;
            FileNamesContainer fileNames;
            fileNames = nameGenerator->GetFileNames( seriesIdentifier );


            reader->SetFileNames( fileNames );

            try
            {
                reader->Update();
            }
            catch (itk::ExceptionObject &ex)
            {
                std::cout << ex << std::endl;
                return EXIT_FAILURE;
            }

            ImageType::Pointer inputImage = reader->GetOutput();


            std::string patientID;
            std::string patientName;
            std::string studyDate;
            std::string aqusitionTime;
            std::string studyDescription;
            std::string seriesDescription;

            std::string manufacturer;
            std::string rescaleIntercept;


            dicomIO->GetValueFromTag("0010|0010", patientName);
            std::cout << patientName << std::endl;

            dicomIO->GetValueFromTag("0010|0020", patientID);
            std::cout << patientID << std::endl;
            
            dicomIO->GetValueFromTag("0008|0020", studyDate);
            std::cout << studyDate << std::endl;

            //dicomIO->GetValueFromTag("0008|0030", studyTime);
            //std::cout << studyTime << std::endl;

            //dicomIO->GetValueFromTag("0008|0031", aqusitionTime);
            //std::cout << aqusitionTime << std::endl;

            dicomIO->GetValueFromTag("0008|0032", aqusitionTime);
            std::cout << aqusitionTime << std::endl;

            dicomIO->GetValueFromTag("0008|1030", studyDescription);
            std::cout << studyDescription << std::endl;

            dicomIO->GetValueFromTag("0008|103e", seriesDescription);
            std::cout << seriesDescription << std::endl;

            dicomIO->GetValueFromTag("0008|0070", manufacturer);
            std::cout << manufacturer << std::endl;

            dicomIO->GetValueFromTag("0028|1052", rescaleIntercept);
            std::cout << rescaleIntercept << std::endl;

            using AddImageFilterType = itk::AddImageFilter <ImageType, ImageType, ImageType>;
            AddImageFilterType::Pointer addImageFilter = AddImageFilterType::New();
            if(rescaleIntercept.length() == 0)
            {
                PixelType newRescaleIntercept = 0;
                refineString(manufacturer);
                if(manufacturer.compare("SIEMENS") == 0)
                {
                    newRescaleIntercept = -1024;
                }
                else
                {
                    newRescaleIntercept = -1000;
                }

                std::cout << "rescale " << newRescaleIntercept << std::endl;

                addImageFilter->SetInput(inputImage);
                addImageFilter->SetConstant2(newRescaleIntercept);
                addImageFilter->Update();

                inputImage = addImageFilter->GetOutput();
            }



            refineString(patientID);
            refineString(studyDate);
            refineString(aqusitionTime);
            refineString(studyDescription);
            refineString(seriesDescription);

            std::string seriesDirectory = argv[2];
            //
            seriesDirectory.append("/");
            seriesDirectory.append(patientID);// + "_" +refineName(patientName));

            itksys::SystemTools::MakeDirectory(seriesDirectory.c_str());
            //
            seriesDirectory.append("/");

            std::string nrrdFilename = seriesDirectory;
            nrrdFilename.append(patientID + "_" + studyDate + "_" + aqusitionTime);
            nrrdFilename.append(".nrrd");

            /*
            if (studyDescription.length() > 0)
            {
                seriesDirectory.append(studyDescription);
                studyDescriptionOld = studyDescription;
            }
            else
            {
                seriesDirectory.append(studyDescriptionOld);
            }
            //
            seriesDirectory.append("/");
            std::string nrrdFilename = seriesDirectory;
            if (seriesDescription.length() > 0)
            {
                seriesDirectory.append(seriesDescription);
                nrrdFilename.append(initialName(patientName) + "_" + seriesDescription + ".nrrd");
            }
            else
            {
                seriesDirectory.append(seriesIdentifier);
                nrrdFilename.append(seriesIdentifier + ".nrrd");
            }
            */




            ////////To write Output images /////////////////////
/*            using WriterType = itk::ImageSeriesWriter< ImageType, DICOMImageType >;

            nameGenerator->SetOutputDirectory( seriesDirectory );
            FileNamesContainer OutputImageNames;
            OutputImageNames = nameGenerator->GetOutputFileNames();

            WriterType::Pointer writer = WriterType::New();
            writer->SetInput(inputImage);
            writer->SetImageIO(dicomIO);
            writer->SetFileNames(OutputImageNames);
            writer->SetMetaDataDictionaryArray(reader->GetMetaDataDictionaryArray());

            std::cout  << "Writing the image as " << std::endl << std::endl;
            std::cout  << seriesDirectory << std::endl << std::endl;*/

            using SingleOutImageWriterType = itk::ImageFileWriter< ImageType>;
            SingleOutImageWriterType::Pointer OutWriter = SingleOutImageWriterType::New();
            OutWriter->SetUseCompression(true);
            OutWriter->SetInput(inputImage);
            OutWriter->SetFileName(nrrdFilename);
            OutWriter->Update();

/*
            try
            {
                writer->Update();
            }
            catch (itk::ExceptionObject &ex)
            {
                std::cout << ex << std::endl;
                return EXIT_FAILURE;
            }
*/
        }
    }
    catch (itk::ExceptionObject &ex)
    {
        std::cout << ex << std::endl;
        return EXIT_FAILURE;
    }

    return EXIT_SUCCESS;
}



<|start_filename|>Tools/ITKUtils/itkGrowCutSegmentationImageFilter.hxx<|end_filename|>
#ifndef itkGrowCutSegmentationImageFilter_txx
#define itkGrowCutSegmentationImageFilter_txx


#include "itkGrowCutSegmentationImageFilter.h"


#include "itkConstNeighborhoodIterator.h"
#include "itkNeighborhoodIterator.h"
#include "itkNeighborhoodAlgorithm.h"
#include "itkImageRegionIterator.h"
#include "itkImageRegionConstIterator.h"
#include "itkConstantBoundaryCondition.h"
#include "itkNumericTraits.h"
#include "itkImageFileWriter.h"
#include "itkLabelMap.h"
#include "itkShapeLabelObject.h"
#include "itkLabelImageToShapeLabelMapFilter.h"

#include "itkProgressReporter.h"
# include "itkIterationReporter.h"

#include <map>
#include <string>
#include <cmath>
#include <vcl_compiler.h>
#include <iostream>
#include <algorithm>
#include <utility>

#include <iostream>
#include <ctime>


// #define LOG_OUTPUT
// #define DEBUG_LEVEL_0

#define USE_LABELSHAPEFILTER

namespace itk
{

template< class TInputImage, class TOutputImage, class TWeightPixelType>
GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>
::GrowCutSegmentationImageFilter()
{
  this->SetNumberOfRequiredInputs(3);

  m_RunOneIteration = false;

  m_SetStateImage = false;

  m_SetDistancesImage = false;

  m_SetMaxSaturationImage = false;

  m_ConfThresh = 0.2;

  m_MaxIterations = 500;
  m_ObjectRadius = 10;

  m_SeedStrength = 1.0;

  m_LabelImage = OutputImageType::New();

  m_WeightImage = WeightImageType::New();

  m_ObjectLabel = static_cast<OutputPixelType>( NumericTraits<OutputPixelType>::max(OutputPixelType() ));
  m_BackgroundLabel = static_cast<OutputPixelType>( 1 );
  m_UnknownLabel = static_cast<OutputPixelType>( NumericTraits<OutputPixelType>::ZeroValue() );

  m_Radius.Fill(1);
}


/**
 * Standard PrintSelf method.
 */
template <class TInputImage, class TOutputImage, class TWeightPixelType>
void
GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>
::PrintSelf(std::ostream& os, Indent indent) const
{
  this->Superclass::PrintSelf(os, indent);
  os << indent << "max iterations : " << m_MaxIterations<< std::endl;
  //   os << indent << "max enemies for attack T1 : " << m_T1<< std::endl;
  // os << indent << "min enemies for submit T2 : " << m_T2<< std::endl;
  os << indent << "starting seed strength :" <<m_SeedStrength<< std::endl;
  //os << indent << "use Algorithm Speed Slow : " << m_UseSlow<< std::endl;
}

template <class TInputImage, class TOutputImage, class TWeightPixelType>
void
GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>
::GenerateInputRequestedRegion()
{
  Superclass::GenerateInputRequestedRegion();
  if ( this->GetInput() )
    {
    InputImagePointer image = const_cast< InputImageType * >( this->GetInput() );
    image->SetRequestedRegionToLargestPossibleRegion();
    }
}

template <class TInputImage, class TOutputImage, class TWeightPixelType>
void GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>
::EnlargeOutputRequestedRegion(DataObject *output)
{
  Superclass::EnlargeOutputRequestedRegion(output);
  output->SetRequestedRegionToLargestPossibleRegion();
}

template <class TInputImage, class TOutputImage, class TWeightPixelType>
void
GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>
::SetStateImage( const OutputImageType *l)
{
  this->ProcessObject::SetNthInput(3, const_cast< OutputImageType *>(l));

  this->SetStateImageOn();
  //m_StateImage->Graft(this->ProcessObject::GetInput(3));
}

template< class TInputImage, class TOutputImage, class TWeightPixelType >
const typename GrowCutSegmentationImageFilter< TInputImage, TOutputImage, TWeightPixelType >::OutputImagePointer
GrowCutSegmentationImageFilter< TInputImage, TOutputImage, TWeightPixelType >
::GetStateImage()
{
  typename OutputImageType::Pointer stateImage = OutputImageType::New();
  stateImage->Graft(static_cast< OutputImageType *>(this->ProcessObject::GetInput(3)));
  return stateImage;
}

template <class TInputImage, class TOutputImage, class TWeightPixelType>
void
GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>
::SetDistancesImage( const WeightImageType* w)
{
  this->ProcessObject::SetNthInput(4, const_cast< WeightImageType*> (w));

  this->SetDistancesImageOn();
  //m_SetDistancesImage = true;
  //m_DistancesImage = static_cast< WeightImageType *>(this->ProcessObject::GetInput(4));
}

template< class TInputImage, class TOutputImage, class TWeightPixelType >
const typename GrowCutSegmentationImageFilter< TInputImage, TOutputImage, TWeightPixelType >::WeightImagePointer
GrowCutSegmentationImageFilter< TInputImage, TOutputImage, TWeightPixelType >
::GetDistancesImage()
{
  typename WeightImageType::Pointer distanceImage = WeightImageType::New();
  distanceImage->Graft(static_cast< WeightImageType *>(this->ProcessObject::GetInput(4)));
  return distanceImage;
  //return const_cast< WeightImageType*> (this->ProcessObject::GetInput(4));
}

template <class TInputImage, class TOutputImage, class TWeightPixelType>
void
GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>
::SetMaxSaturationImage( const WeightImageType* s)
{
  this->ProcessObject::SetNthInput(5, const_cast< WeightImageType *>(s));

  this->SetMaxSaturationImageOn();
  //m_SetMaxSaturationImage = true;
  //m_MaxSaturationImage = static_cast< WeightImageType *>(this->ProcessObject::GetInput(5));
}


template< class TInputImage, class TOutputImage, class TWeightPixelType >
const typename GrowCutSegmentationImageFilter< TInputImage, TOutputImage, TWeightPixelType >::WeightImagePointer
GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>
::GetMaxSaturationImage()
{
  typename WeightImageType::Pointer maxSaturationImage = WeightImageType::New();
  maxSaturationImage->Graft(static_cast< WeightImageType*>(this->ProcessObject::GetInput(5)));
  return maxSaturationImage;
  //return const_cast< WeightImageType*> (this->ProcessObject::GetInput(5));
}

template< class TInputImage, class TOutputImage, class TWeightPixelType >
const typename GrowCutSegmentationImageFilter< TInputImage, TOutputImage, TWeightPixelType >::OutputImagePointer
GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>
::GetLabelImage()
{

  //typename OutputImageType::Pointer labelImage = OutputImageType::New();
  //labelImage->Graft(this->ProcessObject::GetInput(1));
  //return labelImage;
  return m_LabelImage;
  //return const_cast< const OutputImageType*> (this->ProcessObject::GetInput(1));
}

template< class TInputImage, class TOutputImage, class TWeightPixelType >
const typename GrowCutSegmentationImageFilter< TInputImage, TOutputImage, TWeightPixelType >::WeightImagePointer
GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>
::GetStrengthImage()
{
  typename WeightImageType::Pointer weightImage = WeightImageType::New();
  weightImage->Graft(static_cast< WeightImageType*> (this->ProcessObject::GetInput(2)));
  return weightImage;
}


template< class TInputImage, class TOutputImage, class TWeightPixelType >
const typename GrowCutSegmentationImageFilter< TInputImage, TOutputImage, TWeightPixelType >::WeightImagePointer
GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>
::GetUpdatedStrengthImage()
{
  return m_WeightImage;
}

template< class TInputImage, class TOutputImage, class TWeightPixelType >
const typename GrowCutSegmentationImageFilter< TInputImage, TOutputImage, TWeightPixelType >::InputImagePointer
GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>
::GetInputImage()
{
  return const_cast< InputImageType*> (this->GetInput(0));
}

template <class TInputImage, class TOutputImage, class TWeightPixelType>
void
GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>
::Initialize( OutputImageType *output)
{

  // is the output allocated already
  if(output->GetBufferedRegion().GetNumberOfPixels() == 0)
    {
    // allocate memory for the output buffer
    output->SetBufferedRegion( this->GetOutput()->GetRequestedRegion() );
    output->Allocate();
    output->FillBuffer ( m_UnknownLabel );
    }


  if(m_LabelImage->GetBufferedRegion().GetNumberOfPixels() == 0)
    {
    m_LabelImage->CopyInformation( output );
    m_LabelImage->SetBufferedRegion( output->GetBufferedRegion() );
    m_LabelImage->Allocate();
    m_LabelImage->FillBuffer( m_UnknownLabel );
    }


  if( m_WeightImage->GetBufferedRegion().GetNumberOfPixels() == 0)
    {
    m_WeightImage->CopyInformation( output );
    m_WeightImage->SetBufferedRegion( output->GetBufferedRegion() );
    m_WeightImage->Allocate();
    m_WeightImage->FillBuffer( 0.0 );
    }

}

template<class TInputImage, class TOutputImage, class TWeightPixelType>
void GrowCutSegmentationImageFilter< TInputImage, TOutputImage, TWeightPixelType>::
ComputeLabelVolumes(TOutputImage *outputImage, std::vector< unsigned > &volumes,
                    std::vector< unsigned > &physicalVolumes)
{
  //volumes.resize(3);

#ifndef USE_LABELSHAPEFILTER

  std::map<unsigned short, unsigned int>labelMap;
  std::vector< unsigned short > labels;
  unsigned int index = 0;
  ImageRegionConstIteratorWithIndex< OutputImageType > label( outputImage,
                                                              outputImage->GetBufferedRegion());

  for(label.GoToBegin(); !label.IsAtEnd(); ++label)
    {

    unsigned short pix = static_cast<unsigned short>(label.Get());
    std::map<unsigned short, unsigned int>::iterator it = labelMap.find(pix);
    if(it == labelMap.end())
      {
      labelMap.insert(std::pair<unsigned short, unsigned int>(pix, index));
      volumes.push_back(1);
      labels.push_back(pix);
      ++index;
      }
    else
      {
      int i = it->second;
      ++volumes[i];
      }
    }
  // std::cout<<" label \t "<<" volume "<<std::endl;
  // for (unsigned int i = 0; i < index; i++)
  //   {
  //   std::map<unsigned short, unsigned int>::iterator it = labelMap.find(labels[i]);
  //   std::cout<<labels[i]<<"\t"<<volumes[it->second]<<std::endl;
  //   }

#else
  using LabelFilterType = LabelImageToShapeLabelMapFilter< OutputImageType >;
  typename LabelFilterType::Pointer labelFilter = LabelFilterType::New();

  labelFilter->SetInput(outputImage);
  labelFilter->SetBackgroundValue(0);
  labelFilter->Update();

  unsigned long numObjects = labelFilter->GetOutput()->GetNumberOfLabelObjects();
  volumes.resize((const int)numObjects);
  physicalVolumes.resize((const int)numObjects);
  for (unsigned n = 0; n < numObjects; n++)
    {
    volumes[n] = labelFilter->GetOutput()->
      GetNthLabelObject(n)->GetNumberOfPixels();
    physicalVolumes[n] = (unsigned int)
      (labelFilter->GetOutput()->
       GetNthLabelObject(n)->GetNumberOfPixels()/1000.0);
    }
#endif
}


template<class TInputImage, class TOutputImage, class TWeightPixelType>
bool GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>::
InitializeStateImage( OutputImageType *stateImage)
{
  bool setStart = false;
  //  bool setEnd = false;

  typename OutputImageType::Pointer labelImage = OutputImageType::New();
  labelImage->Graft(static_cast< OutputImageType*> (this->ProcessObject::GetInput(1)));

  typename OutputImageType::SizeType lsize =
    labelImage->GetBufferedRegion().GetSize();

  ImageRegionConstIteratorWithIndex< OutputImageType > label( labelImage,
                                                              labelImage->GetBufferedRegion());
  ImageRegionIterator< OutputImageType > state( stateImage,
                                                stateImage->GetBufferedRegion());

  for (unsigned i = 0; i < labelImage->GetImageDimension(); i++){
  m_RoiStart[i] = 0;
  m_RoiEnd[i] = 0;
  }

  bool foundLabels = false;

  for (label.GoToBegin(), state.GoToBegin(); !label.IsAtEnd();
       ++label, ++state)
    {
    if(label.Get() != m_UnknownLabel)
      {
      foundLabels = true;
      OutputIndexType idx = label.GetIndex();

      if(!setStart)
        {
        m_RoiStart = idx;
        m_RoiEnd = idx;
        setStart = true;
        }
      else
        {
        for (unsigned i = 0; i < labelImage->GetImageDimension(); i++)
          {
          if(idx[i] <= m_RoiStart[i])
            {
            m_RoiStart[i] = idx[i];
            }
          if(idx[i] >= m_RoiEnd[i])
            {
            m_RoiEnd[i] = idx[i];
            }
          }
        }

      state.Set(static_cast< OutputPixelType>(LABELED));
      }
    }

  for (unsigned i = 0; i < labelImage->GetImageDimension(); i++)
    {
    /* m_roiStart[i] = (midROI[i]-m_ObjectRadius >=0 && midROI[i]-m_ObjectRadius <= m_roiStart[i]) ?
   (midROI[i]-m_ObjectRadius) : 0; //m_roiStart[i];
   m_roiEnd[i] = (midROI[i]+m_ObjectRadius < (unsigned)lsize[i] && midROI[i]+m_ObjectRadius >= m_roiEnd[i]) ?
   (midROI[i]+m_ObjectRadius) : lsize[i]-1; //m_roiEnd[i];
      */

    m_RoiStart[i] = (m_RoiStart[i] - m_ObjectRadius >= 0)? (m_RoiStart[i] - m_ObjectRadius) : 0;
    m_RoiEnd[i] = (static_cast<unsigned int>(m_RoiEnd[i] + m_ObjectRadius) < lsize[i]) ?
      (m_RoiEnd[i] + m_ObjectRadius) : lsize[i]-1;
    }

  if(!foundLabels)
    return true;

  // std::cout<< " ROI start ";
  // for (unsigned int i = 0; i < labelImage->GetImageDimension(); i++)
  //   {
  //   std::cout<<m_roiStart[i]<<" ";
  //   }
  // std::cout<<" ROI end ";
  // for (unsigned int i = 0; i < labelImage->GetImageDimension(); i++)
  //   {
  //   std::cout<<m_roiEnd[i]<<" ";
  //   }
  // std::cout<<std::endl;

  return false;

}

template<class TInputImage, class TOutputImage, class TWeightPixelType>
void GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>::
InitializeDistancesImage( TInputImage *inputImage,
                          WeightImageType *distances)
{
  using InputNeighborhoodIteratorType = ConstNeighborhoodIterator<InputImageType>;
  using InputRadiusType = typename InputNeighborhoodIteratorType::RadiusType;
  InputRadiusType radiusIn = static_cast<InputRadiusType> (m_Radius);

  InputNeighborhoodIteratorType input(radiusIn,
                                      inputImage,
                                      inputImage->GetBufferedRegion());


  ImageRegionIterator< WeightImageType> dist(distances,
                                             distances->GetBufferedRegion());

  for(input.GoToBegin(), dist.GoToBegin(); !input.IsAtEnd();
      ++input, ++dist)
    {

    WeightPixelType center = static_cast< WeightPixelType>(input.GetCenterPixel());
    WeightPixelType maxDistance = 0.0;
    WeightPixelType distance;

    for (unsigned i = 0; i < input.Size(); i++)
      {
      WeightPixelType pix = static_cast< WeightPixelType>(input.GetPixel(i));

      distance = (pix - center)*(pix-center);
      maxDistance = (distance > maxDistance) ? distance :
        maxDistance;
      }

    dist.Set(maxDistance);
    }
}

template<class TInputImage, class TOutputImage, class TWeightPixelType>
void GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>::
MaskSegmentedImageByWeight( float confThresh)
{
  ImageRegionConstIterator< WeightImageType > weight( m_WeightImage,
                                                      m_WeightImage->GetBufferedRegion());

  ImageRegionIterator< OutputImageType > label( m_LabelImage,
                                                m_LabelImage->GetBufferedRegion());

  //ImageRegionIterator< OutputImageType > label( this->GetOutput(),
  //    this->GetOutput()->GetBufferedRegion());

  weight.GoToBegin();
  label.GoToBegin();
  //  unsigned int numInPixels = 0;
  //unsigned int numOutPixels = 0;

  for(;!weight.IsAtEnd(); ++weight, ++label)
    {
    if(weight.Get() < confThresh)
      {
      label.Set(static_cast< OutputPixelType > (0));
      //  ++numOutPixels;
      }
    //else
    // {
    //  ++numInPixels;
    // }
    }

  //return (numInPixels > numOutPixels);
}


template <class TInputImage, class TOutputImage, class TWeightPixelType>
void
GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>
::GenerateData()
{
  IterationReporter iterate(this, 0, 1);

  // if the filter is configured to run a single iteration, use the superclass
  // GenerateData. This will engage the multi-threaded version of the GenerateData

  typename OutputImageType::Pointer output = this->GetOutput();
  this->Initialize(output);

  if(m_RunOneIteration)
    {
    Superclass::GenerateData();
    iterate.CompletedStep();

    m_LabelImage = this->GetOutput();
    return;
    }


  unsigned int ndims = static_cast< unsigned int>(output->GetImageDimension());

  typename InputImageType::Pointer inputImage = InputImageType::New();
  inputImage->Graft( static_cast< InputImageType*>(this->ProcessObject::GetInput(0)));

  typename OutputImageType::Pointer pixelStateImage = OutputImageType::New();
  typename WeightImageType::Pointer maxDistancesImage = WeightImageType::New();
  typename WeightImageType::Pointer maxSaturationImage = WeightImageType::New();


  bool converged = false;

  // Initialize the state image
  if( !m_SetStateImage)
    {

    pixelStateImage->CopyInformation( this->GetLabelImage() );
    pixelStateImage->SetBufferedRegion( this->GetLabelImage()->GetBufferedRegion() );
    pixelStateImage->Allocate();
    pixelStateImage->FillBuffer(UNLABELED);

    converged = this->InitializeStateImage(pixelStateImage);
    }


  // Initialize the max distance for each pixel in the image
  if( !m_SetDistancesImage)
    {

    maxDistancesImage->CopyInformation(this->GetStrengthImage());
    maxDistancesImage->SetBufferedRegion( this->GetStrengthImage()->GetBufferedRegion() );
    maxDistancesImage->Allocate();
    maxDistancesImage->FillBuffer(UNLABELED);
    this->InitializeDistancesImage(inputImage, maxDistancesImage);
    }


  if( !m_SetMaxSaturationImage )
    {
    maxSaturationImage->CopyInformation( inputImage );
    maxSaturationImage->SetBufferedRegion( inputImage->GetBufferedRegion() );
    maxSaturationImage->Allocate();
    maxSaturationImage->FillBuffer( 0 );
    }

  /////////////////////////////////////////////////////////////////

  // Filter was configured to run until convergence. We need to delegate a different instance of the filter to run on each iteration.
  // std::cout<<" Running filter until convergence .... "<<std::endl;

  typename GrowCutSegmentationImageFilter< TInputImage, TOutputImage, TWeightPixelType >::Pointer
    singleIteration =
    GrowCutSegmentationImageFilter< TInputImage, TOutputImage, TWeightPixelType>::New();
  // set up the singleIteration Filter.
  singleIteration->RunOneIterationOn();
  singleIteration->SetInputImage( this->GetInputImage());
  singleIteration->SetLabelImage( this->GetLabelImage() );
  singleIteration->SetStrengthImage (this->GetStrengthImage() );

  if(m_SetStateImage)
    {
    singleIteration->SetStateImage( this->GetStateImage() );
    }
  else
    {
    singleIteration->SetStateImage(pixelStateImage);
    }

  if(m_SetDistancesImage)
    {
    singleIteration->SetDistancesImage( this->GetDistancesImage() );
    }
  else
    {
    singleIteration->SetDistancesImage(maxDistancesImage);
    }
  if(m_SetMaxSaturationImage)
    {
    singleIteration->SetMaxSaturationImage( this->GetMaxSaturationImage() );
    }
  else
    {
    singleIteration->SetMaxSaturationImage(maxSaturationImage);
    }

  // std::cout<<" Setting other filter parameters ... "<<std::endl;
  singleIteration->SetSeedStrength( this->GetSeedStrength());
  singleIteration->SetMaxIterations( this->GetMaxIterations());
  singleIteration->SetObjectRadius(this->GetObjectRadius() );

  singleIteration->SetROIStart( this->GetROIStart() );
  singleIteration->SetROIEnd( this->GetROIEnd() );

  singleIteration->GetOutput()->SetBufferedRegion(
    this->GetOutput()->GetBufferedRegion() );

  // Create Progress Accumulator for tracking the progress
  //ProgressAccumulator::Pointer progress = ProgressAccumulator::New();
  //progress->SetMiniPipelineFilter( this );
  //progress->RegisterInternalFilter( singleIterationFilter, 1.0f );

  // run until convergence

  unsigned int iter = 0;
  //time_t startTime;
  //time_t endTime;
  time_t start;
  time_t end;

  unsigned int maxRadius = 0;

  unsigned totalROIVolume = 1;

  for (unsigned n = 0; n < ndims; n++)
    {
    //prevUnlabeledPix *= (m_roiEnd[n]-m_roiStart[n]);
    totalROIVolume *= (m_RoiEnd[n]-m_RoiStart[n]);

    if(n == 0)
      {
      maxRadius = m_RoiEnd[n]-m_RoiStart[n];
      }
    else if(maxRadius < (m_RoiEnd[n]-m_RoiStart[n]))
      {
      maxRadius = m_RoiEnd[n]-m_RoiStart[n];
      }
    }

  maxRadius = static_cast< unsigned int>(std::ceil(maxRadius*0.5));

  // float threshSaturation = .96; //.999; // Determine the saturation according to the size of the object
  // float threshUnchanged = 0.05;
  // float threshUnlabeledLimit = 0;

  unsigned threshSaturation = 96; // .96*100
  unsigned threshUnchanged = 5;

  time(&start);

  // set up the grow cut update filter here...

  unsigned prevModifiedPix = 0;
  while (iter < m_MaxIterations && !converged)
    {
    //std::cout<<" running iteration ...... "<<iter<<std::endl;
    singleIteration->Update();
    //std::cout<<" DONE"<<std::endl;
    iterate.CompletedStep();

    // Check for convergence
    //pixelStateImage = singleIteration->GetStateImage();

    unsigned changeablePix =
      singleIteration->GetLabeled();//updateFilter->GetNumLabeledPixels();
    unsigned currLocallySaturatedPix =
      singleIteration->GetLocallySaturated();//updateFilter->GetNumLocallySaturatedPixels();
    unsigned currSaturatedPix =
      singleIteration->GetSaturated();//updateFilter->GetNumSaturatedPixels();

    unsigned unlabeled = (changeablePix + currLocallySaturatedPix + currSaturatedPix <= totalROIVolume) ?
      (static_cast<unsigned > (totalROIVolume) -(currLocallySaturatedPix +
                                                 currSaturatedPix + changeablePix)) : 0;

    unsigned currModified = currLocallySaturatedPix + currSaturatedPix;

    unsigned relUnlabeled = (unlabeled*100)/totalROIVolume;
    unsigned relSaturated = ((currSaturatedPix+currLocallySaturatedPix)*100)/totalROIVolume;

    // unsigned permodified = (currModified > prevModifiedPix) ?
    //  ((currModified-prevModifiedPix)*100)/(std::max(static_cast<unsigned int>(1), std::max(currModified, prevModifiedPix))) :
    //  ((prevModifiedPix - currModified)*100)/(std::max(static_cast<unsigned int>(1), std::max(currModified, prevModifiedPix)));

    converged = (prevModifiedPix == currModified) ||
      (relUnlabeled < threshUnchanged && relSaturated > threshSaturation);

    // if((iter % 10 == 0) || converged)
    //   {
      // std::cout<<" ITER "<<iter<<std::endl;
      // std::cout<<" saturated Pixels "<<currSaturatedPix<<"( "<<
      //   totalROIVolume<<") %"<<currSaturatedPix/(float)totalROIVolume<<std::endl;
      // std::cout<<" locally saturated Pixels "<<currLocallySaturatedPix<<" %"<<
      //   currLocallySaturatedPix/(float)totalROIVolume<<std::endl;
      // std::cout<<"Number of Labeled Pixels "<<changeablePix<<" %"<<
      //   changeablePix/(float)totalROIVolume<<std::endl;
      // std::cout<<"Curr Saturated "<<relSaturated<<" %"<<std::endl;
      // std::cout<<" Prev Modified "<<prevModifiedPix<<" Curr Modified "<<currModified<<std::endl;
      // std::cout<<"Unlabeled pixels "<<unlabeled<<" %"<<
      //   relUnlabeled<<" Per modified "<<permodified<<std::endl;
      // if(converged)
      //   std::cout<<" converged...."<<std::endl;
      // }

    prevModifiedPix = currModified;
    ++iter;
    m_WeightImage = singleIteration->GetUpdatedStrengthImage();

    if(!converged)
      {
      // assign the old output as input
      m_WeightImage->DisconnectPipeline();
      //m_LabelImage->DisconnectPipeline();

      singleIteration->SetStrengthImage(m_WeightImage);
      //singleIteration->SetLabelImage(m_LabelImage);
      singleIteration->GetOutput()->
        SetBufferedRegion(output->GetBufferedRegion() );
      }

    this->UpdateProgress((currSaturatedPix+currLocallySaturatedPix)/(float)totalROIVolume);
    }


  this->UpdateProgress(1.0);
  m_WeightImage = singleIteration->GetUpdatedStrengthImage();
  m_LabelImage = singleIteration->GetOutput();

  this->MaskSegmentedImageByWeight(m_ConfThresh);

  this->GraftOutput(m_LabelImage);

  time(&end);

  std::vector< unsigned > labelVolumes;
  std::vector< unsigned > physicalVolumes;

  this->ComputeLabelVolumes(m_LabelImage, labelVolumes, physicalVolumes);

  // std::cout<<"...................................................."<<std::endl;
  // for (unsigned n = 0; n < labelVolumes.size(); n++)
  //   {
  //   std::cout<<" Object "<<n<<" volume "<<labelVolumes[n]<<" Physical Volume "<<
  //     physicalVolumes[n]<<std::endl;
  //   }

  /*
  using WriterType = ImageFileWriter < OutputImageType >;
  typename WriterType::Pointer writer = WriterType::New();

  std::string outfilename = "segmented_img.mhd";
  writer->SetFileName(outfilename.c_str());
  writer->SetInput(m_LabelImage);
  writer->Update();
  */
  // std::cout<<" Time taken for segmentation "<<iter<<
  //   " : "<<difftime(end, start)<<" secs"<<std::endl;

}


template <class TInputImage, class TOutputImage, class TWeightPixelType>
void GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>
  ::ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread,
                          ThreadIdType itkNotUsed(threadId))
{
  // make a shallow copy of the image to break the pipeline from re-executing in each thread
  typename InputImageType::Pointer inputImage = InputImageType::New();
  inputImage->Graft( static_cast< InputImageType*>(this->ProcessObject::GetInput(0)));

  typename OutputImageType::Pointer labelImage = OutputImageType::New();
  labelImage->Graft(static_cast< OutputImageType*>(this->ProcessObject::GetInput(1)) );

  typename WeightImageType::Pointer weightImage = WeightImageType::New();
  weightImage->Graft(static_cast< WeightImageType*>(this->ProcessObject::GetInput(2)) );

  typename OutputImageType::Pointer stateImage = OutputImageType::New();
  stateImage->Graft(static_cast< OutputImageType*>(this->ProcessObject::GetInput(3)) );

  typename WeightImageType::Pointer distanceImage = WeightImageType::New();
  distanceImage->Graft(static_cast< WeightImageType*>(this->ProcessObject::GetInput(4)) );

  typename WeightImageType::Pointer maxSaturationImage = WeightImageType::New();
  maxSaturationImage->Graft(static_cast< WeightImageType*>(this->ProcessObject::GetInput(5)) );
  //maxSaturationImage->Graft(m_MaxSaturationImage );

  typename OutputImageType::Pointer outputImage = OutputImageType::New();
  outputImage->Graft(this->GetOutput() );


  //ProgressReporter progress(this, threadID, outputRegionForThread.GetNumberOfPixels());

  //////////////////////////////////////////////////////////////////

  typename OutputImageType::IndexType regionStart = outputRegionForThread.GetIndex();
  typename OutputImageType::SizeType regionSize = outputRegionForThread.GetSize();

  //std::cout<<" region "<<regionStart[0]<<" "<<regionStart[1]<<" "<<regionStart[2]<<" "<<threadId<<std::endl;
  //std::cout<<" region size "<<regionSize[0]<<" "<<regionSize[1]<<" "<<regionSize[2]<<" "<<threadId<<std::endl;

  unsigned int Dimension = inputImage->GetImageDimension();

  // Extract the image regions corresponding to the seed and ROI defined by the radius
  typename InputImageType::IndexType iIndex;
  typename InputImageType::SizeType iSize;


  for (unsigned d = 0; d < Dimension; d++)
    {
    unsigned long int x = m_RoiEnd[d]-m_RoiStart[d]+1; // > 0 ? m_roiEnd[d]-m_roiStart[d] : 0;
    iIndex[d] = m_RoiStart[d] > regionStart[d] ? m_RoiStart[d] : regionStart[d];
    iSize[d] = x < regionSize[d] ? (x > 0 ? x : 1) : regionSize[d];

    }

  using InputRegionType = typename InputImageType::RegionType;
  using OutputRegionType = typename OutputImageType::RegionType;
  using WeightRegionType = typename WeightImageType::RegionType;

  InputRegionType input;
  input.SetIndex( iIndex );
  input.SetSize( iSize );

  OutputRegionType label;
  label.SetIndex( iIndex);
  label.SetSize( iSize);

  WeightRegionType weight;
  weight.SetIndex( iIndex);
  weight.SetSize( iSize);

  using InputIteratorType = ConstNeighborhoodIterator<InputImageType>;
  using InputRadiusType = typename InputIteratorType::RadiusType;

  using LabelIteratorType = NeighborhoodIterator<OutputImageType>;
  using LabelRadiusType = typename LabelIteratorType::RadiusType;

  using WeightIteratorType = NeighborhoodIterator<WeightImageType>;
  using WeightRadiusType = typename WeightIteratorType::RadiusType;


  InputRadiusType radiusIn; //= static_cast< InputRadiusType >(m_Radius);
  WeightRadiusType radiusWt; //= static_cast< WeightRadiusType >(m_Radius);
  LabelRadiusType radiusLab; //= static_cast< LabelRadiusType >(m_Radius);
  radiusIn.Fill(1);
  radiusWt.Fill(1);
  radiusLab.Fill(1);

  ConstantBoundaryCondition<InputImageType> cbci;
  cbci.SetConstant( NumericTraits< InputPixelType >::min(InputPixelType()) );

  ConstantBoundaryCondition<OutputImageType> cbcl;
  cbcl.SetConstant( NumericTraits< OutputPixelType >::min(OutputPixelType()) );

  ConstantBoundaryCondition<WeightImageType> cbcw;
  cbcw.SetConstant( NumericTraits< WeightPixelType >::min(WeightPixelType()) );


  InputPixelType minI = NumericTraits< InputPixelType > ::min(InputPixelType());
  WeightPixelType minW = NumericTraits< WeightPixelType > ::min(WeightPixelType());

  /////////////////////////////////////////////////////////////////////////////////////////

  // Boundary Faces Calculator
  using FaceCalculatorType = typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< InputImageType >;
  FaceCalculatorType faceCalculator;
  typename FaceCalculatorType::FaceListType faceList;

  faceList = faceCalculator( inputImage, input, radiusIn);

  // FaceList Iterators
  typename FaceCalculatorType::FaceListType::iterator fitIn;

  //std::cout<<" number faces "<<faceList.size()<<std::endl;
  //bool first = true;
  for( fitIn = faceList.begin(); fitIn != faceList.end(); ++fitIn)
    {

    //fitIn->Print(std::cout);
    //first = true;

    InputIteratorType inputIt; //(radiusIn, inputImage,input);

    LabelIteratorType oldLabelIt; //(radiusLab, labelImage, label);

    using IteratorOut = ImageRegionIterator< OutputImageType >; // newLabelIt(this->GetOutput(), label);

    IteratorOut newLabelIt;

    WeightIteratorType oldWeightIt; //(radiusWt, weightImage, weight);

    using IteratorWeight = ImageRegionIterator< WeightImageType >; //newWeightIt(m_WeightImage, weight);

    IteratorWeight newWeightIt;
    IteratorWeight distancesIt;

    WeightIteratorType maxSatIt;
    LabelIteratorType stateIt;


    inputIt = InputIteratorType(radiusIn, inputImage, *fitIn);

    //oldLabelIt = LabelIteratorType(radiusLab, labelImage, *fitIn);

    oldLabelIt = LabelIteratorType(radiusLab, m_LabelImage, *fitIn);

    oldWeightIt = WeightIteratorType(radiusWt, weightImage, *fitIn);

    maxSatIt = WeightIteratorType(radiusWt, maxSaturationImage, *fitIn);

    stateIt = LabelIteratorType(radiusLab, stateImage, *fitIn);


    inputIt.OverrideBoundaryCondition( &cbci );
    oldLabelIt.OverrideBoundaryCondition( &cbcl );
    stateIt.OverrideBoundaryCondition( &cbcl );
    oldWeightIt.OverrideBoundaryCondition( &cbcw );
    maxSatIt.OverrideBoundaryCondition( &cbcw );

    newLabelIt = IteratorOut(outputImage, *fitIn);
    newWeightIt = IteratorWeight(m_WeightImage, *fitIn);
    distancesIt = IteratorWeight(distanceImage, *fitIn);


    inputIt.GoToBegin(); oldLabelIt.GoToBegin(); newLabelIt.GoToBegin();
    oldWeightIt.GoToBegin(); newWeightIt.GoToBegin(); stateIt.GoToBegin();
    //newstateIt.GoToBegin();
    distancesIt.GoToBegin(); maxSatIt.GoToBegin();

    for (; !inputIt.IsAtEnd(); ++inputIt, ++oldLabelIt,
           ++newLabelIt, ++oldWeightIt, ++newWeightIt,
           ++stateIt, ++distancesIt, ++maxSatIt) // ++newstateIt;
      {

      OutputPixelType s_center = stateIt.GetCenterPixel();
      OutputPixelType l_center = oldLabelIt.GetCenterPixel();
      WeightPixelType w_center = oldWeightIt.GetCenterPixel();

      if(s_center == SATURATED)
        {
        //newstateIt.Set(s_center);
        newLabelIt.Set(l_center);
        newWeightIt.Set(w_center);
        continue;
        }

      InputPixelType f_center = inputIt.GetCenterPixel();

      OutputPixelType winnerLabel = l_center;
      WeightPixelType winnerWeight = w_center;

      WeightPixelType maxDist = distancesIt.Get();

      unsigned int countSaturatedLinks = 0;
      unsigned int countLocalSaturatedLinks = 0;

      bool modified = false;
      WeightPixelType maxWt = 0.0;

      unsigned int nlinks = 0;

      for (unsigned k = 0; k < inputIt.Size(); k++)
        {

        InputPixelType f = inputIt.GetPixel(k);
        WeightPixelType w = oldWeightIt.GetPixel(k);


        if(f == minI && w == minW)
          {
          continue;
          }

        OutputPixelType s = stateIt.GetPixel(k);
        OutputPixelType l = oldLabelIt.GetPixel(k);


        ++nlinks;

        WeightPixelType attackWeight = (f_center - f)*(f_center - f);
        attackWeight = (maxDist > 0) ? (1.0 - attackWeight/maxDist) : 1.0;


        WeightPixelType maxAttackWeight = 0.0; // attackWeight*maxSatIt.GetPixel(k);
        WeightPixelType msat = maxSatIt.GetPixel(k);

        maxAttackWeight = (s == UNLABELED) ? 0.0 :
          ((msat == 0.0) ? (attackWeight * this->GetSeedStrength()) :
           ((s == SATURATED) ? attackWeight * w : attackWeight*msat ) );

        attackWeight *= w;

        maxWt = (maxWt < maxAttackWeight) ? maxAttackWeight : maxWt;


        if(s_center != UNLABELED)
          {
          countSaturatedLinks += (maxAttackWeight <= w_center) ? 1 : 0;
          countLocalSaturatedLinks += (attackWeight <= w_center) ? 1 : 0;
          }

        if(s != UNLABELED && attackWeight > winnerWeight)
          {

          winnerWeight = attackWeight;
          winnerLabel = l;
          modified = true;

          stateIt.SetCenterPixel(LABELED);
          //newstateIt.Set(LABELED);
          }
        }

      if(nlinks > 0)
        {
        if(countSaturatedLinks == nlinks && winnerLabel != m_UnknownLabel)
          {
          stateIt.SetCenterPixel(SATURATED);
          }
        else{
        if (countLocalSaturatedLinks == nlinks && winnerLabel != m_UnknownLabel)
          {
          stateIt.SetCenterPixel(LOCALLY_SATURATED);
          }
        else if (stateIt.GetCenterPixel() != UNLABELED && !modified)
          {
          stateIt.SetCenterPixel(LOCALLY_SATURATED);
          }
        }
        }

      //if(nlinks > 0 && countSaturatedLinks == nlinks && winnerLabel != m_UnknownLabel)
      // {
      //   stateIt.SetCenterPixel(SATURATED);
      //   //newstateIt.Set(SATURATED);
      // }
      //else if((winnerLabel != m_UnknownLabel) &&
      //    ((nlinks > 0) &&
//       (countLocalSaturatedLinks == nlinks)) ||
//      ((stateIt.GetCenterPixel() != UNLABELED) &&
//       !modified))
//      {
//        stateIt.SetCenterPixel(LOCALLY_SATURATED);
//        //newstateIt.Set(LOCALLY_SATURATED);
//      }

      maxSatIt.SetCenterPixel(maxWt);
      newLabelIt.Set(winnerLabel);
      newWeightIt.Set(winnerWeight);
      }
    }

  //std::cout<<" done in thread "<<threadId<<std::endl;
  //time(&end);
  //std::cout<<" time for grow cut update "<<difftime(end,start)<<std::endl;

}


template <class TInputImage, class TOutputImage, class TWeightPixelType>
void
GrowCutSegmentationImageFilter<TInputImage, TOutputImage, TWeightPixelType>
::AfterThreadedGenerateData()
{
  m_Labeled = 0;
  m_LocallySaturated = 0;
  m_Saturated = 0;

  typename OutputImageType::Pointer stateImage = OutputImageType::New();
  stateImage->Graft(static_cast< OutputImageType*>(this->ProcessObject::GetInput(3)) );

  ImageRegionIterator< OutputImageType > state(stateImage, stateImage->GetBufferedRegion());
  ImageRegionIterator< WeightImageType > weight(m_WeightImage, m_WeightImage->GetBufferedRegion());


  for(state.GoToBegin(), weight.GoToBegin(); !state.IsAtEnd(); ++state, ++weight)
    {
    if(state.Get() == LABELED)
      {
      ++m_Labeled;
      }
    if(state.Get() == LOCALLY_SATURATED && weight.Get() >= m_ConfThresh)
      {
      ++m_LocallySaturated;
      }
    if(state.Get() == SATURATED && weight.Get() > m_ConfThresh)
      {
      ++m_Saturated;
      }

    }
}

}//namespace itk

#endif


<|start_filename|>Tools/ITKUtils/itkFastGrowCutSegmentationImageFilter.hxx<|end_filename|>
#ifndef _itkFastGrowCutSegmentationImageFilter_txx_
#define _itkFastGrowCutSegmentationImageFilter_txx_


#include "itkFastGrowCutSegmentationImageFilter.h"


#include "itkConstNeighborhoodIterator.h"
#include "itkNeighborhoodIterator.h"
#include "itkNeighborhoodAlgorithm.h"
#include "itkImageRegionIterator.h"
#include "itkImageRegionConstIterator.h"
#include "itkConstantBoundaryCondition.h"
#include "itkNumericTraits.h"
#include "itkImageFileWriter.h"
#include "itkLabelMap.h"
#include "itkShapeLabelObject.h"
#include "itkLabelImageToShapeLabelMapFilter.h"

#include "itkProgressReporter.h"
# include "itkIterationReporter.h"

#include <vcl_map.h>
#include <vcl_string.h>
#include <vcl_cmath.h>
#include <vcl_algorithm.h>
#include <vcl_utility.h>

#include <iostream>
#include <ctime>

#include "itkTimeProbe.h"

#include "FastGrowCutSegmenter.h"


// #define LOG_OUTPUT
// #define DEBUG_LEVEL_0

#define USE_LABELSHAPEFILTER

namespace itk
{

template< class TInputImage, class TOutputImage>
FastGrowCutSegmentationImageFilter<TInputImage, TOutputImage>
::FastGrowCutSegmentationImageFilter()
{
    this->SetNumberOfRequiredInputs(2);

    m_LabelImage = OutputImageType::New();

    m_fastGC = NULL;
}


/**
 * Standard PrintSelf method.
 */
template <class TInputImage, class TOutputImage>
void
FastGrowCutSegmentationImageFilter<TInputImage, TOutputImage>
::PrintSelf(std::ostream &os, Indent indent) const
{
    this->Superclass::PrintSelf(os, indent);
}

template <class TInputImage, class TOutputImage>
void
FastGrowCutSegmentationImageFilter<TInputImage, TOutputImage>
::GenerateInputRequestedRegion()
{
    Superclass::GenerateInputRequestedRegion();
    if ( this->GetInput() )
    {
        InputImagePointer image = const_cast< InputImageType * >( this->GetInput() );
        image->SetRequestedRegionToLargestPossibleRegion();
    }
}

template <class TInputImage, class TOutputImage>
void FastGrowCutSegmentationImageFilter<TInputImage, TOutputImage>
::EnlargeOutputRequestedRegion(DataObject *output)
{
    Superclass::EnlargeOutputRequestedRegion(output);
    output->SetRequestedRegionToLargestPossibleRegion();
}

template< class TInputImage, class TOutputImage>
const typename FastGrowCutSegmentationImageFilter< TInputImage, TOutputImage>::OutputImagePointer
FastGrowCutSegmentationImageFilter<TInputImage, TOutputImage>
::GetLabelImage()
{

    //typename OutputImageType::Pointer labelImage = OutputImageType::New();
    //labelImage->Graft(this->ProcessObject::GetInput(1));
    //return labelImage;
    return m_LabelImage;
    //return const_cast< const OutputImageType*> (this->ProcessObject::GetInput(1));
}

template< class TInputImage, class TOutputImage>
const typename FastGrowCutSegmentationImageFilter< TInputImage, TOutputImage>::SeedImagePointer
FastGrowCutSegmentationImageFilter<TInputImage, TOutputImage>
::GetSeedImage()
{
    return const_cast< SeedImageType *> (this->GetInput(1));
}

template< class TInputImage, class TOutputImage>
const typename FastGrowCutSegmentationImageFilter< TInputImage, TOutputImage>::InputImagePointer
FastGrowCutSegmentationImageFilter<TInputImage, TOutputImage>
::GetInputImage()
{
    return const_cast< InputImageType *> (this->GetInput(0));
}

template <class TInputImage, class TOutputImage>
void
FastGrowCutSegmentationImageFilter<TInputImage, TOutputImage>
::Initialize( OutputImageType *output)
{

    // is the output allocated already
    if (output->GetBufferedRegion().GetNumberOfPixels() == 0)
    {
        // allocate memory for the output buffer
        output->SetBufferedRegion( this->GetOutput()->GetRequestedRegion() );
        output->Allocate();
    }


    if (m_LabelImage->GetBufferedRegion().GetNumberOfPixels() == 0)
    {
        m_LabelImage->CopyInformation( output );
        m_LabelImage->SetBufferedRegion( output->GetBufferedRegion() );
        m_LabelImage->Allocate();
    }
}

template<class TInputImage, class TOutputImage>
void FastGrowCutSegmentationImageFilter< TInputImage, TOutputImage>::
ComputeLabelVolumes(TOutputImage *outputImage, vcl_vector< unsigned > &volumes,
                    vcl_vector< unsigned > &physicalVolumes)
{
    //volumes.resize(3);

#ifndef USE_LABELSHAPEFILTER

    vcl_map<unsigned short, unsigned int>labelMap;
    vcl_vector< unsigned short > labels;
    unsigned int index = 0;
    ImageRegionConstIteratorWithIndex< OutputImageType > label( outputImage,
            outputImage->GetBufferedRegion());

    for (label.GoToBegin(); !label.IsAtEnd(); ++label)
    {

        unsigned short pix = static_cast<unsigned short>(label.Get());
        vcl_map<unsigned short, unsigned int>::iterator it = labelMap.find(pix);
        if (it == labelMap.end())
        {
            labelMap.insert(vcl_pair<unsigned short, unsigned int>(pix, index));
            volumes.push_back(1);
            labels.push_back(pix);
            ++index;
        }
        else
        {
            int i = it->second;
            ++volumes[i];
        }
    }
    // std::cout<<" label \t "<<" volume "<<std::endl;
    // for (unsigned int i = 0; i < index; i++)
    //   {
    //   vcl_map<unsigned short, unsigned int>::iterator it = labelMap.find(labels[i]);
    //   std::cout<<labels[i]<<"\t"<<volumes[it->second]<<std::endl;
    //   }

#else
    using LabelFilterType = LabelImageToShapeLabelMapFilter< OutputImageType >;
    typename LabelFilterType::Pointer labelFilter = LabelFilterType::New();

    labelFilter->SetInput(outputImage);
    labelFilter->SetBackgroundValue(0);
    labelFilter->Update();

    unsigned long numObjects = labelFilter->GetOutput()->GetNumberOfLabelObjects();
    volumes.resize((const int)numObjects);
    physicalVolumes.resize((const int)numObjects);
    for (unsigned n = 0; n < numObjects; n++)
    {
#if ITK_VERSION_MAJOR < 4
        volumes[n] = labelFilter->GetOutput()->
                     GetNthLabelObject(n)->GetSize();
        physicalVolumes[n] = (unsigned int)
                             (labelFilter->GetOutput()->
                              GetNthLabelObject(n)->GetSize() / 1000.0);
#else
        volumes[n] = labelFilter->GetOutput()->
                     GetNthLabelObject(n)->GetNumberOfPixels();
        physicalVolumes[n] = (unsigned int)
                             (labelFilter->GetOutput()->
                              GetNthLabelObject(n)->GetNumberOfPixels() / 1000.0);
#endif
    }
#endif
}




template <class TInputImage, class TOutputImage>
void
FastGrowCutSegmentationImageFilter<TInputImage, TOutputImage>
::GenerateData()
{
    IterationReporter iterate(this, 0, 1);

    itk::TimeProbe timer;

    m_InitializationFlag = false;
    if (m_fastGC == NULL)
    {
        m_fastGC = new FGC::FastGrowCut<InputPixelType, OutputPixelType>();
    }

    typename OutputImageType::Pointer output = this->GetOutput();
    this->Initialize(output);


    unsigned int ndims = static_cast< unsigned int>(output->GetImageDimension());

    typename InputImageType::Pointer inputImage = InputImageType::New();
    inputImage->Graft( this->ProcessObject::GetInput(0));

    typename SeedImageType::Pointer seedImage = SeedImageType::New();
    seedImage->Graft( this->ProcessObject::GetInput(1));

    timer.Start();

    // Find ROI
    if (!m_InitializationFlag)
    {
        FGC::FindITKImageROI<InputImageType>(inputImage, m_imROI);
        std::cout << "image ROI = [" << m_imROI[0] << "," << m_imROI[1] << "," << m_imROI[2] << ";"  \
                  << m_imROI[3] << "," << m_imROI[4] << "," << m_imROI[5] << "]" << std::endl;
        FGC::ExtractITKImageROI<InputPixelType>(inputImage, m_imROI, m_imSrcVec);
    }

    FGC::ExtractITKImageROI<OutputPixelType>(seedImage, m_imROI, m_imSeedVec);

    // Initialize FastGrowCut
    std::vector<long> imSize(3);
    for (int i = 0; i < 3; i++)
    {
        imSize[i] = m_imROI[i + 3] - m_imROI[i];
    }

    m_fastGC->SetSourceImage(m_imSrcVec);
    m_fastGC->SetSeedlImage(m_imSeedVec);
    m_fastGC->SetImageSize(imSize);
    m_fastGC->SetWorkMode(m_InitializationFlag);

    this->UpdateProgress(0.1);

    // Do Segmentation
    m_fastGC->DoSegmentation();
    //m_fastGC->GetForegroundmage(m_imLabVec);
    m_fastGC->GetLabeImage(m_imLabVec);

    this->UpdateProgress(0.9);

    // Update result
    FGC::UpdateITKImageROI<OutputPixelType>(m_imLabVec, m_imROI, m_LabelImage);

    timer.Stop();

    if (!m_InitializationFlag)
        std::cout << "Initial fast GrowCut segmentation time: " << timer.GetMean() << " seconds\n";
    else
        std::cout << "Adaptive fast GrowCut segmentation time: " << timer.GetMean() << " seconds\n";

    this->UpdateProgress(1.0);

    this->GraftOutput(m_LabelImage);

    vcl_vector< unsigned > labelVolumes;
    vcl_vector< unsigned > physicalVolumes;

    this->ComputeLabelVolumes(m_LabelImage, labelVolumes, physicalVolumes);

    // std::cout<<"...................................................."<<std::endl;
    // for (unsigned n = 0; n < labelVolumes.size(); n++)
    //   {
    //   std::cout<<" Object "<<n<<" volume "<<labelVolumes[n]<<" Physical Volume "<<
    //     physicalVolumes[n]<<std::endl;
    //   }

    /*
    using WriterType = ImageFileWriter < OutputImageType >;
    typename WriterType::Pointer writer = WriterType::New();

    vcl_string outfilename = "segmented_img.mhd";
    writer->SetFileName(outfilename.c_str());
    writer->SetInput(m_LabelImage);
    writer->Update();
    */
    // std::cout<<" Time taken for segmentation "<<iter<<
    //   " : "<<difftime(end, start)<<" secs"<<std::endl;

    if (m_fastGC != NULL)
    {
        delete m_fastGC;
    }

}

}//namespace itk

#endif


<|start_filename|>Tools/GrowCutSegmentation/NoduleSegmentation.cxx<|end_filename|>
#include "itkImage.h"
#include "itkCastImageFilter.h"
#include "itkCommand.h"

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkImageSpatialObject.h"

#include "itkImageRegionConstIterator.h"
#include "itkImageRegionIterator.h"
#include "itkImageRegionIteratorWithIndex.h"

#include "itkBinaryBallStructuringElement.h"
#include "itkBinaryErodeImageFilter.h"
#include "itkBinaryDilateImageFilter.h"
#include "itkBinaryClosingByReconstructionImageFilter.h"
#include "itkBinaryOpeningByReconstructionImageFilter.h"

#include "itkAreaOpeningImageFilter.h"

#include "itkBinaryFillholeImageFilter.h"
#include "itkSliceBySliceImageFilter.h"

#include "itkOtsuMultipleThresholdsImageFilter.h"

#include "itkScalarImageToHistogramGenerator.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkNumericTraits.h"

#include "itkRegionOfInterestImageFilter.h"
#include "itkConnectedComponentImageFilter.h"

#include "itkSpatialObjectToImageFilter.h"
#include "itkEllipseSpatialObject.h"

#include "itkSubtractImageFilter.h"
#include "itkAddImageFilter.h"

#include "itkImageDuplicator.h"

#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkIsotropicResamplerImageFilter.h"

#include "itkSatoVesselnessSigmoidFeatureGenerator.h"

#include "itkGrowCutSegmentationImageFilter.h"

using namespace std;

// Definition of Data type
using InputPixelType = signed short;
using OutputPixelType = unsigned char;
using InternalPixelType = float;


using InImageType = itk::Image<InputPixelType, 3>;
using OutImageType = itk::Image<OutputPixelType, 3>;
using OutImage2DType = itk::Image<OutputPixelType, 2>;
using WeightImageType = itk::Image<InternalPixelType, 3>;


// Definition of IO
using InImageReaderType = itk::ImageFileReader<InImageType>;
using InImageWriterType = itk::ImageFileWriter<InImageType>;
using WeightImageWriterType = itk::ImageFileWriter<WeightImageType>;
using OutImageWriterType = itk::ImageFileWriter<OutImageType>;

using DuplicatorType = itk::ImageDuplicator<OutImageType>;

using RoiFilterType = itk::RegionOfInterestImageFilter<InImageType, InImageType>;


using StructuringElementType = itk::BinaryBallStructuringElement<OutputPixelType, 3>;
using StructuringElement2DType = itk::BinaryBallStructuringElement<OutputPixelType, 2>;
using DilateFilterType = itk::BinaryDilateImageFilter<OutImageType, OutImageType, StructuringElementType>;
using ErodeFilterType = itk::BinaryErodeImageFilter<OutImageType, OutImageType, StructuringElementType>;

using ClosingFilter2DType = itk::BinaryClosingByReconstructionImageFilter<OutImage2DType, StructuringElement2DType>;
using OpeningFilter2DType = itk::BinaryOpeningByReconstructionImageFilter<OutImage2DType, StructuringElement2DType>;


using SliceBySliceFilterType = itk::SliceBySliceImageFilter<OutImageType, OutImageType>;
using FillholeFilter2DType = itk::BinaryFillholeImageFilter<OutImage2DType>;

using AreaOpeningImageFilterType = itk::AreaOpeningImageFilter<OutImageType, OutImageType>;
using ConnectedComponentImageFilterType = itk::ConnectedComponentImageFilter<OutImageType, OutImageType>;

using DiffusionFilterType = itk::CurvatureAnisotropicDiffusionImageFilter<WeightImageType, WeightImageType>;
using IsotropicResamplerType = itk::IsotropicResamplerImageFilter<InImageType, WeightImageType>;
using VesselnessGeneratorType = itk::SatoVesselnessSigmoidFeatureGenerator<3>;
using SpatialObjectType = VesselnessGeneratorType::SpatialObjectType;
using InputImageSpatialObjectType = VesselnessGeneratorType::InputImageSpatialObjectType;
using OutputImageSpatialObjectType = itk::ImageSpatialObject<VesselnessGeneratorType::Dimension, InternalPixelType>;


using ThresholdFilterType = itk::BinaryThresholdImageFilter<WeightImageType, OutImageType>;
using ConnectedComponentImageFilterType = itk::ConnectedComponentImageFilter<OutImageType, OutImageType>;

using GrowCutFilterType = itk::GrowCutSegmentationImageFilter<WeightImageType, OutImageType>;

using SubFilterType = itk::SubtractImageFilter<OutImageType, OutImageType, OutImageType>;

// To print the progress
class ShowProgressObject
{
public:
    ShowProgressObject(itk::ProcessObject * o)
    {
        m_Process = o;
    }

    void
    ShowProgress()
    {
        cout << "\rProgress: "
             << static_cast<unsigned int>(100.0 * m_Process->GetProgress()) << "%";
        if (m_Process->GetProgress() == 1) {
            cout << endl;
        }
    }

    itk::ProcessObject::Pointer m_Process;
};


int
main(int argc, char * argv[])
{
    if (argc < 8) {
        cerr << "Missing Parameters " << endl;
        cerr << "Usage: " << argv[0];
        cerr << " InputImageFile SeedPoint[x y z] NoduleSize[long short] OutputImageFile" << endl;
        return EXIT_FAILURE;
    }

    const unsigned int x  = atoi(argv[2]);
    const unsigned int y  = atoi(argv[3]);
    const int sliceNumber = atoi(argv[4]);
    double longD  = atof(argv[5]);
    double shortD = atof(argv[6]);

    double objectSize = longD;
    double priorSegmentStrength = 0.003;
    double contrastNoiseRatio   = 0.8;

    double shortR = shortD / 2;
    double longR  = longD / 2;

    unsigned char foreground = 255;
    unsigned char background = 0;


    /////////To Read CT images ////////////////////
    InImageReaderType::Pointer InImageReader = InImageReaderType::New();
    InImageReader->SetFileName(argv[1]);
    InImageReader->Update();


    InImageType::Pointer image = InImageReader->GetOutput();


    if (contrastNoiseRatio > 1.0) {
        contrastNoiseRatio /= 100.0;
    }

    if (priorSegmentStrength > 1.0) {
        priorSegmentStrength /= 100.0;
    }

    // set ROI
    InImageType::IndexType index     = image->GetLargestPossibleRegion().GetIndex();
    InImageType::SizeType size       = image->GetLargestPossibleRegion().GetSize();
    InImageType::SpacingType spacing = image->GetSpacing();
    InImageType::SizeType roiSize;
    InImageType::PointType center_phy;

    roiSize[0] = static_cast<unsigned int>((longD + 10) / spacing[0]);
    roiSize[1] = static_cast<unsigned int>((longD + 10) / spacing[1]);
    roiSize[2] = static_cast<unsigned int>((longD + 10) / spacing[2]);

    OutImageType::IndexType roiStart;
    OutImageType::IndexType roiEnd;

    OutImageType::IndexType center;

    center[0] = x;
    center[1] = y;

    if (sliceNumber > 0)
        center[2] = size[2] - sliceNumber - 1;
    else
        center[2] = -sliceNumber;

    cout << center << endl;
    image->TransformIndexToPhysicalPoint(center, center_phy);
    cout << center_phy << endl;


    roiStart[0] = center[0] - roiSize[0] / 2;
    roiStart[1] = center[1] - roiSize[1] / 2;
    roiStart[2] = center[2] - roiSize[2] / 2;
    roiEnd[0]   = roiStart[0] + roiSize[0];
    roiEnd[1]   = roiStart[1] + roiSize[1];
    roiEnd[2]   = roiStart[2] + roiSize[2];

    cout << roiSize << endl;
    cout << roiStart << endl;
    cout << roiEnd << endl;


    unsigned int ndims = image->GetImageDimension();

    cout << " objectSize (radius) " << objectSize << endl;
    OutImageType::PixelType radius = static_cast<OutImageType::PixelType>(objectSize);

    /////////////////////////////////////////////////////////////////

    InImageType::IndexType istart;
    InImageType::SizeType isize;

    for (unsigned n = 0; n < ndims; n++) {
        roiStart[n] = roiStart[n] < 0 ? 0 : roiStart[n];
        istart[n]   = roiStart[n];
        isize[n]    = roiEnd[n] - roiStart[n];
        isize[n]    = (static_cast<unsigned int>(isize[n] + istart[n]) <= size[n]) ?
          isize[n] : (size[n] - istart[n] - 1);
    }

    cout << " roiStart " << roiStart << " roiEnd " << roiEnd << " " << endl;
    cout << " istart " << istart << " isize " << isize << " " << endl;


    InImageType::RegionType iRegion;
    iRegion.SetIndex(istart);
    iRegion.SetSize(isize);

    RoiFilterType::Pointer fInput = RoiFilterType::New();
    fInput->SetRegionOfInterest(iRegion);
    fInput->SetInput(image);
    fInput->Update();

    InImageType::Pointer inImage = InImageType::New();
    inImage = fInput->GetOutput();


    {
        InImageWriterType::Pointer writer = InImageWriterType::New();
        writer->SetUseCompression(true);
        writer->SetInput(inImage);
        writer->SetFileName(string(argv[7]) + "in.nrrd");
        writer->Update();
    }


    /////////////////////////////////////////////////////////////////

    WeightImageType::SpacingType outSpacing = image->GetSpacing();
    WeightImageType::Pointer inImageIso;
    if (outSpacing[2] != outSpacing[0]) {
		if(outSpacing[2] > outSpacing[0])
			outSpacing[2] = outSpacing[0];
		else {
			outSpacing[0] = outSpacing[2];
			outSpacing[1] = outSpacing[2];
		}

        cout << "Out spacing " << outSpacing << endl;

        IsotropicResamplerType::Pointer isotropicResampler = IsotropicResamplerType::New();

        isotropicResampler->SetInput(inImage);
        isotropicResampler->SetOutputSpacing(outSpacing);
        isotropicResampler->Update();

        inImageIso = isotropicResampler->GetOutput();
    } else {
        using CastImageFilterType = itk::CastImageFilter<InImageType, WeightImageType>;
        CastImageFilterType::Pointer castImageFilter = CastImageFilterType::New();
        castImageFilter->SetInput(inImage);
        castImageFilter->Update();

        inImageIso = castImageFilter->GetOutput();
    }
    iRegion = inImageIso->GetLargestPossibleRegion();


    {
        WeightImageWriterType::Pointer writer = WeightImageWriterType::New();

        writer->SetUseCompression(true);
        writer->SetInput(inImageIso);
        writer->SetFileName(string(argv[7]) + "in-iso.nrrd");
        writer->Update();
    }

    // Definition of Morphological operation sturcture elements
    StructuringElementType::RadiusType elementRadius;

    StructuringElement2DType structuringElement2D1;
    structuringElement2D1.SetRadius(1);
    structuringElement2D1.CreateStructuringElement();

    StructuringElementType structuringElement1;
    structuringElement1.SetRadius(1);
    structuringElement1.CreateStructuringElement();

    /*elementRadius.Fill(static_cast<unsigned int>(round(3 / outSpacing[0])));

    StructuringElementType structuringElement3;
    structuringElement3.SetRadius(elementRadius);
    structuringElement3.CreateStructuringElement();


	elementRadius.Fill(static_cast<unsigned int>(round(5 / outSpacing[0])));

	StructuringElementType structuringElement5;
	structuringElement5.SetRadius(elementRadius);
	structuringElement5.CreateStructuringElement();*/

    elementRadius.Fill(static_cast<unsigned int>(round(10 / outSpacing[0])));

    StructuringElementType structuringElement10;
    structuringElement10.SetRadius(elementRadius);
    structuringElement10.CreateStructuringElement();

    elementRadius.Fill(static_cast<unsigned int>(round(15 / outSpacing[0])));

    StructuringElementType structuringElement15;
    structuringElement15.SetRadius(elementRadius);
    structuringElement15.CreateStructuringElement();


    OutImageType::IndexType roiCenter;
    inImageIso->TransformPhysicalPointToIndex(center_phy, roiCenter);
    for (unsigned n = 0; n < ndims; n++) {
        roiStart[n] = iRegion.GetIndex(n);
        roiEnd[n]   = roiStart[n] + iRegion.GetSize(n) - 1;
    }

    cout << " roiStart " << roiStart << " roiEnd " << roiEnd << " " << endl;
    cout << " roiCenter " << roiCenter << endl;


    OutImageType::Pointer labImage   = OutImageType::New();
    WeightImageType::Pointer wtImage = WeightImageType::New();

    labImage->CopyInformation(inImageIso);
    labImage->SetBufferedRegion(inImageIso->GetBufferedRegion());
    labImage->Allocate();
    labImage->FillBuffer(0);

    wtImage->CopyInformation(inImageIso);
    wtImage->SetBufferedRegion(inImageIso->GetBufferedRegion());
    wtImage->Allocate();
    wtImage->FillBuffer(0);


	DiffusionFilterType::Pointer diffusionFilter = DiffusionFilterType::New();
	// Diffusion filter parameters
	diffusionFilter->SetInput(inImageIso);
	diffusionFilter->SetNumberOfIterations(5);
	diffusionFilter->SetTimeStep(outSpacing[0] / 16);
	diffusionFilter->SetConductanceParameter(2);
	diffusionFilter->UseImageSpacingOn();
	diffusionFilter->Update();

    {
		WeightImageWriterType::Pointer writer = WeightImageWriterType::New();
		writer->SetUseCompression(true);
		writer->SetInput(diffusionFilter->GetOutput());
		writer->SetFileName(string(argv[7]) + "filt.nrrd");
		writer->Update();
    }
	inImageIso = diffusionFilter->GetOutput();


    DilateFilterType::Pointer dilateFilter10 = DilateFilterType::New();
    DilateFilterType::Pointer dilateFilter15 = DilateFilterType::New();

    /////////////////////////////////////////////////////////////////
    {
        InImageType::IndexType start;
        InImageType::SizeType size;
        InImageType::RegionType centerRegion;

        size.Fill(round(max(shortD/ outSpacing[0], 5.0)));
        start[0] = roiCenter[0] - round(max(shortD / 2 / outSpacing[0],2.0));
        start[1] = roiCenter[1] - round(max(shortD / 2 / outSpacing[1],2.0));
        start[2] = roiCenter[2] - round(max(shortD / 2 / outSpacing[2],2.0));
        centerRegion.SetSize(size);
        centerRegion.SetIndex(start);

        cout << centerRegion << endl;

        // itk::ImageRegionIterator< InImageType > centerPixels(image, centerRegion);
        // itk::ImageRegionIterator< InImageType > centerPixels(inImage, centerRegion);
        // itk::ImageRegionIterator< WeightImageType > centerPixels(diffusionFilter->GetOutput(), centerRegion);
        itk::ImageRegionIterator<WeightImageType> centerPixels(inImageIso, centerRegion);
		const double th_low = -800;
		const double th_high = 200;
        double Tc        = -300;
        double stdTc     = 250;
		double pTc = 0;
		double pStdTc = 0;
		double varTc = 0;
        double numPixels = 125;
        int iter         = 0;
		int numMaxPix = 0;

        while (	numPixels > 30 &&
				numMaxPix < 100 &&
				iter++ < 20 && 
				stdTc > 200 && 
				(pTc != Tc || pStdTc != stdTc)) {
            double sumPix  = 0;
            double sumPix2 = 0;
            double Tc2     = 0;

            numPixels = 0;
			numMaxPix = 0;

            centerPixels.GoToBegin();
            while (!centerPixels.IsAtEnd()) {
                double pixel = centerPixels.Get();
                if (pixel > th_high){
					numMaxPix++;
				}
                if ((pixel > th_low  && pixel < th_high) && (pixel > Tc - 1 * stdTc && pixel < Tc + 2 * stdTc)) {
                    sumPix  = sumPix + pixel;
                    sumPix2 = sumPix2 + pixel * pixel;
					numPixels++;
                    //cout << pixel << " ";
                } else {
                    //cout << "(" << pixel << ") ";
                }

                ++centerPixels;
            }
            pTc    = Tc;
            pStdTc = stdTc;
            if(numPixels>10){
              Tc     = sumPix / numPixels;
              Tc2    = sumPix2 / numPixels;
              varTc  = Tc2 - Tc * Tc;
              stdTc  = sqrt(varTc);
            }
            if(stdTc == 0) stdTc = 50;

            if (Tc < th_low || Tc > th_high || stdTc > 400  || stdTc < 50) {
                if (Tc > th_high) Tc = th_high-250;
				if (Tc < th_low) Tc = th_low+250;
                stdTc = 100;
                cout << endl << "Reset initial: Tc=" << Tc << " stdTc=" << stdTc << endl;
                continue;
            }

            cout << endl;
            cout << "--------" << endl;
            cout << numPixels << "     => " << Tc << "+-" << stdTc << endl;
        }
        cout << "Threshold: " << Tc - 2 * stdTc << endl;

        ThresholdFilterType::Pointer thresholdFilter    = ThresholdFilterType::New();
        ThresholdFilterType::Pointer thresholdFilter_bg = ThresholdFilterType::New();
        // thresholdFilter->SetInput(inImage);
        // thresholdFilter_bg->SetInput(inImage);
        thresholdFilter->SetInput(inImageIso);
        thresholdFilter_bg->SetInput(inImageIso);

        thresholdFilter->SetLowerThreshold(Tc - 2 * stdTc);
        thresholdFilter->SetUpperThreshold(5000);

        thresholdFilter_bg->SetLowerThreshold(th_low);
        thresholdFilter_bg->SetUpperThreshold(5000);

        // thresholdFilter_bg->SetLowerThreshold(Tc - 200);
        // thresholdFilter_bg->SetUpperThreshold(Tc + 400);

        thresholdFilter->SetInsideValue(1);
        thresholdFilter->SetOutsideValue(0);
        thresholdFilter->Update();

        thresholdFilter_bg->SetInsideValue(0);
        thresholdFilter_bg->SetOutsideValue(1);
        thresholdFilter_bg->Update();

        OutImageType::Pointer thresholdImage = thresholdFilter->GetOutput();
		{
            SliceBySliceFilterType::Pointer sliceBySliceFilter = SliceBySliceFilterType::New();

            FillholeFilter2DType::Pointer fillholeFilter2D = FillholeFilter2DType::New();
            fillholeFilter2D->SetFullyConnected(true);
            fillholeFilter2D->SetForegroundValue(1);

            sliceBySliceFilter->SetInput(thresholdImage);
            sliceBySliceFilter->SetFilter(fillholeFilter2D);

            sliceBySliceFilter->Update();

            thresholdImage = sliceBySliceFilter->GetOutput();
        }

        /*{
            SliceBySliceFilterType::Pointer sliceBySliceFilter = SliceBySliceFilterType::New();

            OpeningFilter2DType::Pointer openingFilter2D = OpeningFilter2DType::New();
            openingFilter2D->SetKernel(structuringElement2D1);
            openingFilter2D->SetForegroundValue(1);

            sliceBySliceFilter->SetInput(thresholdImage);
            sliceBySliceFilter->SetFilter(openingFilter2D);

            sliceBySliceFilter->Update();

            thresholdImage = sliceBySliceFilter->GetOutput();
        }


        {
            AreaOpeningImageFilterType::Pointer openingFilter = AreaOpeningImageFilterType::New();
            openingFilter->SetInput(thresholdImage);
            openingFilter->SetLambda(3);
            // openingFilter->FullyConnectedOn();
            openingFilter->Update();
            cout << "Opening" << endl;
            thresholdImage = openingFilter->GetOutput();
        }*/
        cout << "Thresholding" << endl;

        // Fixed seed
		{
        	InImageType::IndexType start;
	        InImageType::SizeType size;
	        InImageType::RegionType centerRegion;

        	size.Fill(3);
	        start[0] = roiCenter[0] - 1;
	        start[1] = roiCenter[1] - 1;
	        start[2] = roiCenter[2] - 1;
	        centerRegion.SetSize(size);
	        centerRegion.SetIndex(start);
	        
            itk::ImageRegionIterator<OutImageType> centerPixels(thresholdImage, centerRegion);
			
            centerPixels.GoToBegin();
            while (!centerPixels.IsAtEnd()) {
				centerPixels.Set(1);

                ++centerPixels;
            }
        }

		/////////////////////////////////////////////////////////////////
		cout << "Connected Component" << endl;
		double target = 0;
		ConnectedComponentImageFilterType::Pointer connected = ConnectedComponentImageFilterType::New();
		try
		{
			AreaOpeningImageFilterType::Pointer openingFilter = AreaOpeningImageFilterType::New(); // delete small particles
			openingFilter->SetInput(thresholdImage);
			openingFilter->SetLambda(3);
			openingFilter->FullyConnectedOn();
			openingFilter->Update();
			cout << "Opening" << endl;
			
			connected->SetInput(openingFilter->GetOutput());
			connected->Update();


			/////////////////////////////////////////////////////////////////
			cout << "Select Center Component" << endl;
			unsigned int count[255] = { 0, };
			{
				itk::ImageRegionIterator<OutImageType> centerPixels(connected->GetOutput(), centerRegion);
				itk::ImageRegionIterator<OutImageType> centerPixels1(thresholdFilter_bg->GetOutput(), centerRegion);

				centerPixels.GoToBegin();
				centerPixels1.GoToBegin();
				while (!centerPixels.IsAtEnd()) {
					unsigned char label = centerPixels.Get();
					if (label > 0) {
						++count[label];
						centerPixels1.Set(0);
					}

					++centerPixels;
					++centerPixels1;
				}
			}
			unsigned int max_count = 0;
			for (unsigned c = 1; c < 255; c++) {
				if (count[c] > max_count) {
					max_count = count[c];
					target = c;
					cout << c << ", " << max_count << endl;
				}
			}
			cout << target << endl;

		}
		catch (itk::ExceptionObject &excp)
		{
			cerr << "Exception thrown while writing the series " << endl;
			cerr << excp << endl;
		}

        OutImageType::Pointer targetImage = OutImageType::New();
        int numSeeds = 0;

        targetImage->CopyInformation(inImageIso);
        targetImage->SetBufferedRegion(inImageIso->GetBufferedRegion());
        targetImage->Allocate();
        targetImage->FillBuffer(0);
        {
            itk::ImageRegionIterator<OutImageType> label(connected->GetOutput(), centerRegion);
            itk::ImageRegionIterator<OutImageType> targetLabel(targetImage, centerRegion);
            label.GoToBegin();
            targetLabel.GoToBegin();
            while (!label.IsAtEnd()) {
                if (label.Get() == target) {
                    targetLabel.Set(1);
                    ++numSeeds;
                } else {
					targetLabel.Set(0); 
				}
                ++label;
                ++targetLabel;
            }
        }

        cout << numSeeds << endl;



        /*{
            SubFilterType::Pointer subFilter = SubFilterType::New();
            subFilter->SetInput1(targetImage);
            subFilter->SetInput2(thresholdFilter_bg->GetOutput());
            subFilter->Update();

            targetImage = subFilter->GetOutput();

            itk::ImageRegionIterator<OutImageType> targetLabel(targetImage, targetImage->GetBufferedRegion());
            targetLabel.GoToBegin();
            while (!targetLabel.IsAtEnd()) {
                if (targetLabel.Get() != 1) targetLabel.Set(0);
                ++targetLabel;
            }

            AreaOpeningImageFilterType::Pointer openingFilter = AreaOpeningImageFilterType::New();
            openingFilter->SetInput(targetImage);
            openingFilter->SetLambda(5);
            openingFilter->Update();
            targetImage = openingFilter->GetOutput();
        }*/


        ErodeFilterType::Pointer erodeFilterB = ErodeFilterType::New();
        erodeFilterB->SetInput(thresholdFilter_bg->GetOutput());
        erodeFilterB->SetKernel(structuringElement1);
        erodeFilterB->SetErodeValue(1);
        erodeFilterB->Update();


        dilateFilter10->SetInput(targetImage);
        dilateFilter10->SetKernel(structuringElement10);
        dilateFilter10->SetDilateValue(1);
        dilateFilter10->Update();

        dilateFilter15->SetInput(targetImage);
        dilateFilter15->SetKernel(structuringElement15);
        dilateFilter15->SetDilateValue(1);
        dilateFilter15->Update();
        {
            ConnectedComponentImageFilterType::Pointer connected = ConnectedComponentImageFilterType::New();
            connected->SetInput(dilateFilter15->GetOutput());
            connected->Update();

            OutImageType::Pointer image = connected->GetOutput();

            itk::ImageRegionIterator<OutImageType> centerLabels(image, centerRegion);

            DuplicatorType::Pointer duplicator = DuplicatorType::New();
            duplicator->SetInputImage(image);
            duplicator->Update();
            OutImageType::Pointer tempImage = duplicator->GetOutput();


            double target = connected->GetOutput()->GetPixel(roiCenter);

            OutImageType::SizeType roiSize;
            roiSize[0] = static_cast<unsigned int>(longD / outSpacing[0]);
            roiSize[1] = static_cast<unsigned int>(longD / outSpacing[1]);
            roiSize[2] = static_cast<unsigned int>(longD / outSpacing[2]);

            OutImageType::IndexType roiStart;
            OutImageType::IndexType roiEnd;

            roiStart[0] = roiCenter[0] - roiSize[0] / 2;
            roiStart[1] = roiCenter[1] - roiSize[1] / 2;
            roiStart[2] = roiCenter[2] - roiSize[2] / 2;
            roiEnd[0]   = roiStart[0] + roiSize[0];
            roiEnd[1]   = roiStart[1] + roiSize[1];
            roiEnd[2]   = roiStart[2] + roiSize[2];

            cout << roiSize << endl;
            cout << roiStart << endl;
            cout << roiEnd << endl;


            itk::ImageRegionIterator<OutImageType> tempLabel(tempImage, tempImage->GetBufferedRegion());
            itk::ImageRegionIteratorWithIndex<OutImageType> label(dilateFilter15->GetOutput(),
              image->GetBufferedRegion());
            label.GoToBegin();
            tempLabel.GoToBegin();
            while (!label.IsAtEnd()) {
                OutImageType::IndexType idx = label.GetIndex();

                if (idx[0] >= roiStart[0] && idx[1] >= roiStart[1] && idx[2] >= roiStart[2] &&
                  idx[0] < roiEnd[0] && idx[1] < roiEnd[1] && idx[2] < roiEnd[2])
                {
                    if (tempLabel.Get() == target) label.Set(1);
                } else { 
					label.Set(0); 
				}
                ++label;
                ++tempLabel;
            }
        }


        // generate seed areas
        {
            // itk::ImageRegionIterator< OutImageType > bg(erodeFilterB->GetOutput(), erodeFilterB->GetOutput()->GetBufferedRegion());
            itk::ImageRegionIterator<OutImageType> bg(thresholdFilter_bg->GetOutput(),
              thresholdFilter_bg->GetOutput()->GetBufferedRegion());
            itk::ImageRegionIterator<OutImageType> gray(dilateFilter10->GetOutput(),
              dilateFilter10->GetOutput()->GetBufferedRegion());
            itk::ImageRegionIterator<OutImageType> truebg(dilateFilter15->GetOutput(),
              dilateFilter15->GetOutput()->GetBufferedRegion());

            itk::ImageRegionIterator<OutImageType> targetLabel(targetImage, targetImage->GetBufferedRegion());
            itk::ImageRegionIterator<OutImageType> newLabel(labImage, labImage->GetBufferedRegion());
            itk::ImageRegionIterator<WeightImageType> weight(wtImage, wtImage->GetBufferedRegion());

            targetLabel.GoToBegin();
            newLabel.GoToBegin();
            weight.GoToBegin();

            bg.GoToBegin();
            gray.GoToBegin();
            truebg.GoToBegin();
            while (!targetLabel.IsAtEnd()) {
                weight.Set(0.0);

                if (targetLabel.Get() == 1) { // taget
                    weight.Set(contrastNoiseRatio);
                    newLabel.Set(1);
                }

				if (bg.Get() == 1) { // background
                    weight.Set(contrastNoiseRatio);
                    newLabel.Set(2);
                }
				else if (gray.Get() == 0) {
					weight.Set(contrastNoiseRatio);
					newLabel.Set(3);
				}

                if (truebg.Get() == 0) { // definetely background, opposite value (0 is bg)
                    //weight.Set(contrastNoiseRatio);
                    weight.Set(priorSegmentStrength);
                    newLabel.Set(4);
                }

                ++targetLabel;
                ++newLabel;
                ++weight;
                ++bg;
                ++gray;
                ++truebg;
            }
        }

        {
            OutImageWriterType::Pointer writer = OutImageWriterType::New();
            writer->SetUseCompression(true);
            writer->SetInput(labImage);
            writer->SetFileName(string(argv[7]) + "Tc-label.nrrd");
            writer->Update();
        }
    }

    cout << "label image generation completed" << endl;


    /////////////////////////////////////////////////////////////////

    // GrowCut Segmentation
    GrowCutFilterType::Pointer filter = GrowCutFilterType::New();

    // Segmentation filter parameters
    filter->SetInput(inImageIso);
    filter->SetLabelImage(labImage);

    filter->SetStrengthImage(wtImage);

    filter->SetSeedStrength(contrastNoiseRatio);
    filter->SetObjectRadius((unsigned int) objectSize);

    //
    ShowProgressObject progressWatch(filter);
    itk::SimpleMemberCommand<ShowProgressObject>::Pointer command;
    command = itk::SimpleMemberCommand<ShowProgressObject>::New();
    command->SetCallbackFunction(&progressWatch,
      &ShowProgressObject::ShowProgress);
    filter->AddObserver(itk::ProgressEvent(), command);

    filter->Update();


    OutImageType::Pointer outputImageROI = OutImageType::New();

    outputImageROI = filter->GetOutput();

    cout << "Done running segmentation filter " << endl;

    /////////////////////////////////////////////////////////////////

    // Post processing
    {
        itk::ImageRegionIterator<OutImageType> label(outputImageROI, outputImageROI->GetBufferedRegion());
        itk::ImageRegionIterator<OutImageType> truebg(dilateFilter15->GetOutput(), dilateFilter15->GetOutput()->GetBufferedRegion());

        label.GoToBegin();
        truebg.GoToBegin();
        while (!label.IsAtEnd()) {
            if (label.Get() != 1) label.Set(0);
            if (truebg.Get() == 0) label.Set(0);
            ++label;
            ++truebg;
        }

        SliceBySliceFilterType::Pointer sliceBySliceFilter = SliceBySliceFilterType::New();

        OpeningFilter2DType::Pointer openingFilter2D = OpeningFilter2DType::New();
        openingFilter2D->SetKernel(structuringElement2D1);
        openingFilter2D->SetForegroundValue(1);

        sliceBySliceFilter->SetInput(outputImageROI);
        sliceBySliceFilter->SetFilter(openingFilter2D);

        sliceBySliceFilter->Update();

        outputImageROI = sliceBySliceFilter->GetOutput();

        AreaOpeningImageFilterType::Pointer openingFilter = AreaOpeningImageFilterType::New();
        openingFilter->SetInput(outputImageROI);
        openingFilter->SetLambda(3);
        openingFilter->FullyConnectedOn();
        openingFilter->Update();
        cout << "Opening" << endl;
        outputImageROI = openingFilter->GetOutput();
    }

    /*{
     *  OutImageWriterType::Pointer writer = OutImageWriterType::New();
     *  writer->SetUseCompression(true);
     *  writer->SetInput(outputImageROI);
     *  writer->SetFileName(string(argv[7]) + "mask-label.nrrd");
     *  writer->Update();
     * }*/


    /////////////////////////////////////////////////////////////////

    /*OutImageType::Pointer outputImage = OutImageType::New();
     * // allocate outputImage first
     * outputImage->CopyInformation(image);
     * outputImage->SetSpacing(spacing);
     * outputImage->SetBufferedRegion(image->GetBufferedRegion());
     * outputImage->Allocate();
     * outputImage->FillBuffer(0);
     *
     * itk::ImageRegionIterator< OutImageType > filterOut(outputImageROI, outputImageROI->GetBufferedRegion());
     * itk::ImageRegionIterator< OutImageType > out(outputImage, oRegion);
     *
     * for (filterOut.GoToBegin(), out.GoToBegin(); !filterOut.IsAtEnd(); ++filterOut, ++out)
     * {
     *  if (filterOut.Get() == 1)
     *      out.Set(1);
     *  //out.Set(filterOut.Get());
     * }
     *
     * cout << "writeImage_spacing = " << outputImage->GetSpacing() << endl;
     * cout << "writeImage_origin  = " << outputImage->GetOrigin() << endl;
     * cout << "writeImage_LargestPossibleRegion = " << outputImage->GetLargestPossibleRegion() << endl;*/


    ////////To write Output images /////////////////////
    OutImageWriterType::Pointer OutImageWriter = OutImageWriterType::New();
    OutImageWriter->SetUseCompression(true);
    OutImageWriter->SetInput(outputImageROI);
    OutImageWriter->SetFileName(argv[7]);

    try
    {
        OutImageWriter->Update();
    }
    catch (itk::ExceptionObject &excp)
    {
        cerr << "Exception thrown while writing the series " << endl;
        cerr << excp << endl;
        return EXIT_FAILURE;
    }


    cout << "Saved the tumor mask image" << endl;

    // Vesselness feature generate
    {
        VesselnessGeneratorType::Pointer vesselnessFeatureGenerator = VesselnessGeneratorType::New();
        InputImageSpatialObjectType::Pointer inputSpatialObject     = InputImageSpatialObjectType::New();

        using CastImageFilterType = itk::CastImageFilter<WeightImageType, InImageType>;
        CastImageFilterType::Pointer castImageFilter = CastImageFilterType::New();
        castImageFilter->SetInput(inImageIso);

        inputSpatialObject->SetImage(castImageFilter->GetOutput());

        vesselnessFeatureGenerator->SetInput(inputSpatialObject);
        vesselnessFeatureGenerator->SetSigma(1.0);
        vesselnessFeatureGenerator->SetAlpha1(0.1);
        vesselnessFeatureGenerator->SetAlpha2(2.0);
        vesselnessFeatureGenerator->SetSigmoidAlpha(-10.0);
        vesselnessFeatureGenerator->SetSigmoidBeta(40.0);
        vesselnessFeatureGenerator->SetUseVesselEnhancingDiffusion(true);

        vesselnessFeatureGenerator->Update();

        cout << "vesselness feature generation completed" << endl;

        typename SpatialObjectType::Pointer vesselness =
          const_cast<SpatialObjectType *>(vesselnessFeatureGenerator->GetFeature());
        typename OutputImageSpatialObjectType::Pointer outputObject =
          dynamic_cast<OutputImageSpatialObjectType *>( vesselness.GetPointer() );
        typename WeightImageType::Pointer vesselImage =
          const_cast<WeightImageType *>(outputObject->GetImage());
        vesselImage->DisconnectPipeline();

    itk::ImageRegionIterator<WeightImageType> vesselnessit(vesselImage, vesselImage->GetBufferedRegion());
    vesselnessit.GoToBegin();
    while (!vesselnessit.IsAtEnd()) {
      vesselnessit.Set(1 - vesselnessit.Get());

      ++vesselnessit;
    }

        {
            WeightImageWriterType::Pointer writer = WeightImageWriterType::New();
            writer->SetUseCompression(true);
            writer->SetInput(vesselImage);
            writer->SetFileName(string(argv[7]) + "-vessel.nrrd");
            writer->Update();
        }
    }

    return EXIT_SUCCESS;
} // main


<|start_filename|>Tools/ITKUtils/itkMaskedSpatialObjectToImageFilter.h<|end_filename|>
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkMaskedSpatialObjectToImageFilter.h,v $
  Language:  C++
  Date:      $Date: 2006/07/19 15:21:11 $
  Version:   $Revision: 1.00 $

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

     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.

=========================================================================*/
#ifndef __itkMaskedSpatialObjectToImageFilter_h
#define __itkMaskedSpatialObjectToImageFilter_h

#include "itkImageSource.h"
#include "itkSpatialObjectToImageFilter.h"
#include "itkConceptChecking.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkNumericTraits.h"
#include "itkImageRegionConstIterator.h"
#include "itkProgressReporter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkBinaryDilateImageFilter.h"
#include "itkBinaryBallStructuringElement.h"
#include "itkResampleImageFilter.h"
#include "itkCenteredAffineTransform.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkBinaryThresholdImageFilter.h"

namespace itk
{

/** \class MaskedSpatialObjectToImageFilter
 * \brief Base class for filters that take a SpatialObject
 *        as input and produce an image as output, using a Mask to
 *        help decrease computation time.
 *  The user can set a pre-defined mask using SetMask(), else this
 *  filter will automatically compute a mask at a lower resolution.
 *
 *  See "Using a Mask to Decrease Computation Time for SpatialObject
 *  to Image Conversions", Insight Journal, 2006 July-December.
 */
template <class TInputSpatialObject, class TOutputImage>
class ITK_EXPORT MaskedSpatialObjectToImageFilter : public SpatialObjectToImageFilter<TInputSpatialObject, TOutputImage>
{
public:
  /** Standard class type alias. */
  using Self = MaskedSpatialObjectToImageFilter;
  using Superclass = SpatialObjectToImageFilter<TInputSpatialObject, TOutputImage>;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;
  using SizeType = typename TOutputImage::SizeType;
  using PointType = typename TOutputImage::PointType;
  using OutputImageType = TOutputImage;
  using OutputPixelType = typename OutputImageType::PixelType;
  using OutputImagePointer = typename OutputImageType::Pointer;
  using ValueType = typename OutputImageType::ValueType;
  using SpacingType = typename OutputImageType::SpacingType;

  /** Method for creation through the object factory. */
  itkNewMacro(Self);

  /** Run-time type information (and related methods). */
  itkTypeMacro(MaskedSpatialObjectToImageFilter, SpatialObjectToImageFilter);

  /** Superclass type alias. */
  using OutputImageRegionType = typename Superclass::OutputImageRegionType;

  /** Some type alias for input SpatialObject. */
  using InputSpatialObjectType = TInputSpatialObject;
  using InputSpatialObjectPointer = typename InputSpatialObjectType::Pointer;
  using InputSpatialObjectConstPointer = typename InputSpatialObjectType::ConstPointer;
  using ChildrenListType = typename TInputSpatialObject::ChildrenListType;

  /** Some convenient type alias. */
  using MaskCoordRep = double;
  using MaskPixelType = unsigned char;
  using MaskImageType = itk::Image<MaskPixelType, itkGetStaticConstMacro(OutputImageDimension)>;
  using MaskImagePointer = typename MaskImageType::Pointer;
  using MaskImageConstPointer = typename MaskImageType::ConstPointer;

  /** Set/Get the image input of this process object.  */
  virtual void SetInput(const InputSpatialObjectType *object);
  const InputSpatialObjectType *GetInput(void);

  /** Set/Get the input mask image to help speed up processing.
    * If a mask is not specified, one will be automatically generated
    * using GetMaskResampleFactor(). */
  virtual void SetMask(const MaskImageType *mask);
  const MaskImageType *GetMask(void);

  /** If the Mask image is not set, this filter will automatically
   *  generate a mask by creating low-resolution version of the SpatialObject.
   *  This parameter sets the resampling factor for all dimensions.
   *  Example: 4.0 means that the mask is generated at 1/4 the output size.
   *  A large the factor means faster mask generation and less accuracy.
   *  A smaller factor means slower mask generation and more accuarcy.
   *  Default = 4.0 */
  itkSetMacro(MaskResampleFactor, double);
  itkGetMacro(MaskResampleFactor, double);

  /** Sets/Gets the size of the binary dilation structuring element used to
   *  generate the Mask.
   *  The larger the value, the more accurate the mask but the slower the generation.
   *  Default = 2. */
  itkSetMacro(MaskDilationSize, unsigned int);
  itkGetMacro(MaskDilationSize, unsigned int);

protected:
  MaskedSpatialObjectToImageFilter();
  ~MaskedSpatialObjectToImageFilter();

  virtual void GenerateData();
  virtual void GenerateDataUsingMask();
  virtual void SetupOutputImage();
  virtual void GenerateMask();

  bool m_UserSuppliedMask;
  double m_MaskResampleFactor;
  unsigned int m_MaskDilationSize;

  virtual void PrintSelf(std::ostream &os, Indent indent) const;

private:
  MaskedSpatialObjectToImageFilter(const Self &); //purposely not implemented
  void operator=(const Self &);                   //purposely not implemented
};

} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkMaskedSpatialObjectToImageFilter.hxx"
#endif

#endif


<|start_filename|>Tools/ContourTools/ExtractBoundary.cxx<|end_filename|>
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkResampleImageFilter.h>
#include <itkCastImageFilter.h>
#include <itkCommand.h>

#include <itkRegionOfInterestImageFilter.h>
#include <itkExtractImageFilter.h>

// Operations on Images

#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIterator.h>
#include <vnl/vnl_math.h>

// For Feature
#include <itkRescaleIntensityImageFilter.h>


#include <itkMinimumMaximumImageCalculator.h>
#include <itkStatisticsImageFilter.h>

#include <itkScalarImageToTextureFeaturesFilter.h>
#include <itkScalarImageToRunLengthFeaturesFilter.h>

#include <itkLabelGeometryImageFilter.h>
#include <itkConnectedComponentImageFilter.h>
#include <itkRelabelComponentImageFilter.h>
#include <itkLabelStatisticsImageFilter.h>

// For Label Object Representation
#include <itkLabelObject.h>
#include <itkLabelMap.h>
#include <itkShapeLabelObject.h>
#include <itkStatisticsLabelObject.h>
#include <itkBinaryImageToLabelMapFilter.h>
#include <itkLabelMapToLabelImageFilter.h>
#include <itkLabelImageToLabelMapFilter.h>
#include <itkShapeLabelMapFilter.h>
#include <itkStatisticsLabelMapFilter.h>
#include <itkBinaryImageToShapeLabelMapFilter.h>
#include <itkBinaryImageToStatisticsLabelMapFilter.h>

// For island removing
#include <itkBinaryShapeKeepNObjectsImageFilter.h>

#include <itkBinaryFillholeImageFilter.h>
#include <itkSliceBySliceImageFilter.h>
#include <itkSubtractImageFilter.h>
#include <itkAndImageFilter.h>

#include "itkFlatStructuringElement.h"
#include "itkBinaryErodeImageFilter.h"
#include "itkBinaryDilateImageFilter.h"

// For file operation
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>

// For threshold
#include <itkThresholdImageFilter.h>
#include <itkBinaryThresholdImageFilter.h>

#include "ITKUtils.h"


// Definition of Feature Extraction
using MaxMinFilterType = itk::MinimumMaximumImageCalculator< InputImageType >;

using namespace std;

int main( int argc, char *argv[] )
{


    if ( argc < 3)
    {
        cerr << "Missing Parameters " << endl;
        cerr << "Usage = " << argv[0];
        cerr << " InputLabelImage OutputLabelImage ThicknessOfBoundary={2} Mode={2D}" << endl;
        return EXIT_FAILURE;
    }


    double radius = 2;
    if (argc > 3)
    {
        radius = atoi(argv[3]);
    }

    char mode3D = 0;
    if (argc > 4)
    {
        cout << argv[4] << endl;
        if(strcmp(argv[4],"3D")==0)
            mode3D = 1;
        else if(strcmp(argv[4],"3D1")==0)
            mode3D = 2;
    }

    bool isTumorOnly = 0;
    bool isExpand = 0;
    MaskPixelType tumorLabel = 307;
    if (argc > 5)
    {
        cout << argv[5] << endl;
        if(strcmp(argv[5],"T")==0)
            isTumorOnly = 1;
        if(strcmp(argv[5],"T1")==0)
        {
            isTumorOnly = 1;
            tumorLabel = 1;
        }
        if(strcmp(argv[5],"TE")==0)
        {
            isTumorOnly = 1;
            isExpand = 1;
        }
    }


    string inputLabelImageName = argv[1];
    string outputLabelImageName = argv[2];


    cout << "Input Label Image Name = " << inputLabelImageName << endl;
    cout << "Output Label Image Name = " << outputLabelImageName << endl;


    ///////// To Read the Mask image ////////////////////
    cout << "Get the mask from Input" << endl;
    LabelImageType::Pointer inputLabelImage = ReadImageFile<LabelImageType>(inputLabelImageName);

    // To check coordinate of the input label image
    SpacingType inputImageSpacing = inputLabelImage->GetSpacing();
    OriginType  inputImageOrigin  = inputLabelImage->GetOrigin();
    RegionType  inputImageRegion  = inputLabelImage->GetLargestPossibleRegion();
    SizeType    inputImageSize    = inputImageRegion.GetSize();


    /////////////////////////////

    cout << "Input Image Spacing = " << inputImageSpacing << endl;
    cout << "Input Image Origin = " << inputImageOrigin << endl;
    cout << "Input Image Size = " << inputImageSize << endl << endl << endl;


    // Image ROI
    MaskImageType::RegionType maskRegion;

    MaskImageType::Pointer normalMaskImage;
    MaskImageType::Pointer tumorMaskImage;

    normalMaskImage = GetMaskImage(inputLabelImage, 306); // normal
    tumorMaskImage = GetMaskImage(inputLabelImage, tumorLabel); // tumor
    maskRegion = GetRoi(tumorMaskImage);
    ExpandRoi(tumorMaskImage, maskRegion, radius+2);
    BoundingCheck(normalMaskImage, inputLabelImage, maskRegion, inputImageRegion);
    BoundingCheck(tumorMaskImage, inputLabelImage, maskRegion, inputImageRegion);

    MaskImageType::IndexType roiStart;
    MaskImageType::IndexType roiEnd;

    RegionToIndex(maskRegion, roiStart, roiEnd);

    MaskImageType::Pointer tumorMaskImageROI = ApplyRoi<MaskImageType>(tumorMaskImage, maskRegion);
    MaskImageType::Pointer normalMaskImageROI = ApplyRoi<MaskImageType>(normalMaskImage, maskRegion);


    //Definition of Morphological operation sturcture elements
    using StructuringElementType = itk::FlatStructuringElement<Dimension>;
    using DilateFilterType = itk::BinaryDilateImageFilter<MaskImageType, MaskImageType, StructuringElementType>;
    using ErodeFilterType = itk::BinaryErodeImageFilter<MaskImageType, MaskImageType, StructuringElementType>;

    using AndImageFilterType = itk::AndImageFilter <MaskImageType>;
    using SubtractImageFilterType = itk::SubtractImageFilter <MaskImageType, MaskImageType>;

    using MaskImage2DType = itk::Image<MaskPixelType, Dimension-1>;
    using StructuringElement2DType = itk::FlatStructuringElement<Dimension-1>;
    using DilateFilter2DType = itk::BinaryDilateImageFilter<MaskImage2DType, MaskImage2DType, StructuringElement2DType>;
    using ErodeFilter2DType = itk::BinaryErodeImageFilter<MaskImage2DType, MaskImage2DType, StructuringElement2DType>;

    using SliceBySliceFilterType = itk::SliceBySliceImageFilter<MaskImageType, MaskImageType>;

    StructuringElementType::RadiusType elementRadius;
    elementRadius.Fill(static_cast<unsigned int>(round(radius/inputImageSpacing[0])));
    elementRadius[2] = static_cast<unsigned int>(round(radius/inputImageSpacing[2]));
    if(mode3D == 2)
        elementRadius[2] = 1;

    StructuringElement2DType::RadiusType elementRadius2D;
    elementRadius2D.Fill(static_cast<unsigned int>(round(radius/inputImageSpacing[0])));

    StructuringElementType structuringElement = StructuringElementType::Ball(elementRadius, true);
    StructuringElement2DType structuringElement2D = StructuringElement2DType::Ball(elementRadius2D, true);


    MaskImageType::Pointer outputNormalMaskImageROI;
    MaskImageType::Pointer outputTumorMaskImageROI;
    {
        // Nb = N&dilate(T,size), 306
        AndImageFilterType::Pointer andFilter = AndImageFilterType::New();
        itk::ImageToImageFilter<MaskImageType, MaskImageType>::Pointer filterPointer;
        if(mode3D==0)
        {
            cout << elementRadius2D << endl;
            SliceBySliceFilterType::Pointer sliceBySliceFilter = SliceBySliceFilterType::New();

            DilateFilter2DType::Pointer dilate2DFilter = DilateFilter2DType::New();
            dilate2DFilter->SetKernel(structuringElement2D);
            dilate2DFilter->SetDilateValue(maskValue);

            sliceBySliceFilter->SetInput(tumorMaskImageROI);
            sliceBySliceFilter->SetFilter(dilate2DFilter);

            filterPointer = sliceBySliceFilter;
        }
        else
        {
            cout << elementRadius << endl;
            DilateFilterType::Pointer dilateFilter = DilateFilterType::New();

            dilateFilter->SetInput(tumorMaskImageROI);
            dilateFilter->SetKernel(structuringElement);
            dilateFilter->SetDilateValue(maskValue);

            filterPointer = dilateFilter;
        }
        ShowProgressObject progressWatch(filterPointer);
        itk::SimpleMemberCommand<ShowProgressObject>::Pointer command;
        command = itk::SimpleMemberCommand<ShowProgressObject>::New();
        command->SetCallbackFunction(&progressWatch,
                                     &ShowProgressObject::ShowProgress);

        filterPointer->AddObserver(itk::ProgressEvent(), command);


        filterPointer->Update();

        if(isTumorOnly)
        {
            outputNormalMaskImageROI = MaskImageType::New();
            outputNormalMaskImageROI->CopyInformation( normalMaskImageROI );
            outputNormalMaskImageROI->SetSpacing(inputImageSpacing);
            outputNormalMaskImageROI->SetBufferedRegion( normalMaskImageROI->GetBufferedRegion() );
            outputNormalMaskImageROI->Allocate();
            outputNormalMaskImageROI->FillBuffer(0);
            outputTumorMaskImageROI = filterPointer->GetOutput();
        }
        else
        {
            andFilter->SetInput1(normalMaskImageROI);
            andFilter->SetInput2(filterPointer->GetOutput());
            andFilter->Update();

            outputNormalMaskImageROI = andFilter->GetOutput();
        }
        cout << endl;
    }

    if(isTumorOnly == 0)
    {
        // Tb = T & dilate(N,size), 307
        AndImageFilterType::Pointer andFilter = AndImageFilterType::New();
        itk::ImageToImageFilter<MaskImageType, MaskImageType>::Pointer filterPointer;
        if(mode3D==0)
        {
            cout << elementRadius2D << endl;
            SliceBySliceFilterType::Pointer sliceBySliceFilter = SliceBySliceFilterType::New();

            DilateFilter2DType::Pointer dilate2DFilter = DilateFilter2DType::New();
            dilate2DFilter->SetKernel(structuringElement2D);
            dilate2DFilter->SetDilateValue(maskValue);

            sliceBySliceFilter->SetInput(normalMaskImageROI);
            sliceBySliceFilter->SetFilter(dilate2DFilter);

            filterPointer = sliceBySliceFilter;
        }
        else
        {
            DilateFilterType::Pointer dilateFilter = DilateFilterType::New();

            dilateFilter->SetInput(normalMaskImageROI);
            dilateFilter->SetKernel(structuringElement);
            dilateFilter->SetDilateValue(maskValue);

            filterPointer = dilateFilter;
        }
        ShowProgressObject progressWatch(filterPointer);
        itk::SimpleMemberCommand<ShowProgressObject>::Pointer command;
        command = itk::SimpleMemberCommand<ShowProgressObject>::New();
        command->SetCallbackFunction(&progressWatch,
                                     &ShowProgressObject::ShowProgress);

        filterPointer->AddObserver(itk::ProgressEvent(), command);


        filterPointer->Update();

        if(isTumorOnly==0) // intersection
        {
            andFilter->SetInput1(tumorMaskImageROI);
            andFilter->SetInput2(filterPointer->GetOutput());
            andFilter->Update();
            outputTumorMaskImageROI = andFilter->GetOutput();
        }

        cout << endl;
    }

    cout << "Complete Extraction!" << endl;

    LabelImageType::Pointer outputLabelImage = LabelImageType::New();
    // allocate outputLabelImage first
    outputLabelImage->CopyInformation( inputLabelImage );
    outputLabelImage->SetSpacing(inputImageSpacing);
    outputLabelImage->SetBufferedRegion( inputLabelImage->GetBufferedRegion() );
    outputLabelImage->Allocate();
    outputLabelImage->FillBuffer(0);


    if(isTumorOnly == 0)
    {
        itk::ImageRegionIterator< MaskImageType > normalOut(outputNormalMaskImageROI, outputNormalMaskImageROI->GetBufferedRegion());
        itk::ImageRegionIterator< MaskImageType > tumorOut(outputTumorMaskImageROI, outputTumorMaskImageROI->GetBufferedRegion());
        itk::ImageRegionIterator< LabelImageType > out(outputLabelImage, inputImageRegion);

        for (normalOut.GoToBegin(), tumorOut.GoToBegin(), out.GoToBegin(); !normalOut.IsAtEnd(); ++normalOut, ++tumorOut, ++out)
        {
            if(normalOut.Get() > 0 )
            {
                //cout << "N ";
                out.Set(306);
            }
            if(tumorOut.Get() > 0 )
            {
                //cout << "T ";
                out.Set(307);
            }
        }
    }
    else
    {
        itk::ImageRegionIterator< MaskImageType > tumor(tumorMaskImageROI, tumorMaskImageROI->GetBufferedRegion());
        itk::ImageRegionIterator< MaskImageType > tumorOut(outputTumorMaskImageROI, outputTumorMaskImageROI->GetBufferedRegion());
        itk::ImageRegionIterator< LabelImageType > out(outputLabelImage, inputImageRegion);

        for (tumor.GoToBegin(), tumorOut.GoToBegin(), out.GoToBegin(); !tumor.IsAtEnd(); ++tumor, ++tumorOut, ++out)
        {
            if((isExpand == 0 && tumor.Get() ==0 && tumorOut.Get() > 0)
                || (isExpand==1 && tumorOut.Get() > 0))
            {
                //cout << "T ";
                out.Set(tumorLabel);
            }
        }
    }




    std::cout << "writeImage_spacing = " << outputLabelImage->GetSpacing() << std::endl;
    std::cout << "writeImage_origin  = "  << outputLabelImage->GetOrigin() << std::endl;
    std::cout << "writeImage_LargestPossibleRegion = " << outputLabelImage->GetLargestPossibleRegion() << std::endl;


    ////////To write Output images /////////////////////
    try
    {
        WriteImageFile<LabelImageType>(outputLabelImage,outputLabelImageName);
    }
    catch (itk::ExceptionObject &excp)
    {
        std::cerr << "Exception thrown while writing the series " << std::endl;
        std::cerr << excp << std::endl;
        return EXIT_FAILURE;
    }


    std::cout << "Saved the normal tissue and tumor boundary image" << std::endl;

    return EXIT_SUCCESS;
}


<|start_filename|>Tools/DICOMTools/DICOM2NRRDConverter.cxx<|end_filename|>
#include "itkImage.h"
#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkImageFileReader.h"
#include "itkImageSeriesReader.h"
#include "itkImageSeriesWriter.h"

#include "itkImageFileWriter.h"
#include "itkTextOutput.h"

#include "itksys/SystemTools.hxx"

#include <boost/algorithm/string.hpp>


using ShortPixelType = signed short;
using UShortPixelType = unsigned short;
using FloatPixelType = float;
using ShortImageType = itk::Image< ShortPixelType, 3 >;
using ShortDICOMImageType = itk::Image< ShortPixelType, 2 >;
using UShortImageType = itk::Image< UShortPixelType, 3 >;
using UShortDICOMImageType = itk::Image< UShortPixelType, 2 >;
using FloatImageType = itk::Image< FloatPixelType, 3 >;
using FloatDICOMImageType = itk::Image< FloatPixelType, 2 >;

using NamesGeneratorType = itk::GDCMSeriesFileNames;
using FileNamesContainer = std::vector< std::string >;

int refineString(std::string &in)
{
	boost::trim(in);
	boost::replace_all(in, "/", "_");
	boost::replace_all(in, " ", "_");

	return in.length();
}

std::string refineName(std::string &fullName)
{
	std::vector<std::string> strs;
	std::string refineName;

	boost::algorithm::to_lower(fullName);
	boost::split(strs, fullName, boost::is_any_of("^"));

	boost::trim_if(strs[0], !boost::algorithm::is_alpha());
	boost::trim_if(strs[1], !boost::algorithm::is_alpha());

	strs[0][0] = toupper(strs[0][0]);
	strs[1][0] = toupper(strs[1][0]);

	refineName = strs[1] + "_" + strs[0];

	return refineName;
}

std::string initialName(std::string &fullName)
{
	std::vector<std::string> strs;
	std::string initialName;

	boost::algorithm::to_lower(fullName);
	boost::split(strs, fullName, boost::is_any_of("^"));

	boost::trim_if(strs[0], !boost::algorithm::is_alpha());
	boost::trim_if(strs[1], !boost::algorithm::is_alpha());

	strs[0][0] = toupper(strs[0][0]);
	strs[1][0] = toupper(strs[1][0]);

	initialName = strs[1][0] + strs[0];

	return initialName;
}

std::string getModality(std::string fileName)
{
	using ReaderType = itk::ImageFileReader<ShortDICOMImageType>;
	using ImageIOType = itk::GDCMImageIO;

	ImageIOType::Pointer dicomIO = ImageIOType::New();
	ReaderType::Pointer reader = ReaderType::New();
	reader->SetImageIO(dicomIO);

	reader->SetFileName(fileName);
	reader->Update();

	std::string modality;

	dicomIO->GetValueFromTag("0008|0060", modality);
	//std::cout << modality << std::endl;

	return modality;
}

template <typename TImageType>
int readAndWrite(FileNamesContainer fileNames, std::string seriesDirectory)
{
	using ReaderType = itk::ImageSeriesReader<TImageType>;
	using ImageIOType = itk::GDCMImageIO;

	ImageIOType::Pointer dicomIO = ImageIOType::New();
	typename ReaderType::Pointer reader = ReaderType::New();
	reader->SetImageIO(dicomIO);


	reader->SetFileNames(fileNames);

	try
	{
		reader->Update();
	}
	catch (itk::ExceptionObject &ex)
	{
		std::cout << ex << std::endl;
		return EXIT_FAILURE;
	}


	std::string patientID;
	std::string patientName;
	std::string studyDate;
	std::string studyTime;
	std::string seriesTime;
	std::string aqusitionTime;
	std::string modality;
	std::string studyDescription;
	std::string seriesDescription;
	std::string timepointID;
	std::string convolutionKernel;

	dicomIO->GetValueFromTag("0010|0010", patientName);
	std::cout << patientName << std::endl;

	dicomIO->GetValueFromTag("0010|0020", patientID);
	std::cout << patientID << std::endl;


	dicomIO->GetValueFromTag("0008|0020", studyDate);
	std::cout << studyDate << std::endl;

	dicomIO->GetValueFromTag("0008|0030", studyTime);
	std::cout << studyTime << std::endl;

	dicomIO->GetValueFromTag("0008|0031", seriesTime);
	std::cout << seriesTime << std::endl;

	dicomIO->GetValueFromTag("0008|0032", aqusitionTime);
	std::cout << aqusitionTime << std::endl;

	dicomIO->GetValueFromTag("0008|0060", modality);
	std::cout << modality << std::endl;

	dicomIO->GetValueFromTag("0008|1030", studyDescription);
	std::cout << studyDescription << std::endl;

	dicomIO->GetValueFromTag("0008|103e", seriesDescription);
	std::cout << seriesDescription << std::endl;

	dicomIO->GetValueFromTag("0012|0050", timepointID);
	std::cout << timepointID << std::endl;

	dicomIO->GetValueFromTag("0018|1210", convolutionKernel);
	std::cout << convolutionKernel << std::endl;

	refineString(patientID);
	refineString(studyDate);
	refineString(aqusitionTime);
	refineString(studyDescription);
	refineString(seriesDescription);
	refineString(convolutionKernel);


	//
	seriesDirectory.append("/");
	seriesDirectory.append(patientID);// + "_" +refineName(patientName));

	itksys::SystemTools::MakeDirectory(seriesDirectory.c_str());
	//
	seriesDirectory.append("/");

	std::string nrrdFilename = seriesDirectory;
	nrrdFilename.append(patientID + "_" + studyDate + "_" + modality);
	/*if (seriesDescription.length() > 0)
	{
		nrrdFilename.append(seriesDescription);
	}
	else
	{
		if (convolutionKernel.length() > 0)
			nrrdFilename.append(convolutionKernel);
		else
			nrrdFilename.append(studyTime);
	}*/

	nrrdFilename.append(".nrrd");

	/*
	if (studyDescription.length() > 0)
	{
	seriesDirectory.append(studyDescription);
	studyDescriptionOld = studyDescription;
	}
	else
	{
	seriesDirectory.append(studyDescriptionOld);
	}
	//
	seriesDirectory.append("/");
	std::string nrrdFilename = seriesDirectory;
	if (seriesDescription.length() > 0)
	{
	seriesDirectory.append(seriesDescription);
	nrrdFilename.append(initialName(patientName) + "_" + seriesDescription + ".nrrd");
	}
	else
	{
	seriesDirectory.append(seriesIdentifier);
	nrrdFilename.append(seriesIdentifier + ".nrrd");
	}
	*/




	////////To write Output images /////////////////////
	//using WriterType = itk::ImageSeriesWriter<ImageType>;

	//nameGenerator->SetOutputDirectory(seriesDirectory);
	//FileNamesContainer OutputImageNames;
	//OutputImageNames = nameGenerator->GetOutputFileNames();

	//WriterType::Pointer writer = WriterType::New();
	//writer->SetInput(reader->GetOutput());
	//writer->SetImageIO(dicomIO);
	//writer->SetFileNames(OutputImageNames);
	//writer->SetMetaDataDictionaryArray(reader->GetMetaDataDictionaryArray());

	std::cout << "Writing the image as " << std::endl << std::endl;
	std::cout << seriesDirectory << std::endl << std::endl;

	using SingleOutImageWriterType = itk::ImageFileWriter<TImageType>;
	typename SingleOutImageWriterType::Pointer OutWriter = SingleOutImageWriterType::New();
	OutWriter->SetUseCompression(true);
	OutWriter->SetInput(reader->GetOutput());
	OutWriter->SetFileName(nrrdFilename);
	OutWriter->Update();

	return 0;
}

int main(int argc, char *argv[])
{
	if (argc < 3)
	{
		std::cerr << "Usage: " << std::endl;
		std::cerr << argv[0] << " DicomDirectory OutputDicomDirectory"
			<< std::endl;
		return EXIT_FAILURE;
	}

	itk::OutputWindow::SetInstance(itk::TextOutput::New());

	NamesGeneratorType::Pointer nameGenerator = NamesGeneratorType::New();
	nameGenerator->SetUseSeriesDetails(true);
	nameGenerator->AddSeriesRestriction("0008|0021");
	nameGenerator->SetRecursive(true);
	nameGenerator->SetDirectory(argv[1]);

	std::string seriesDirectory = argv[2];

	try
	{
		using SeriesIdContainer = std::vector< std::string >;
		const SeriesIdContainer &seriesUID = nameGenerator->GetSeriesUIDs();

		if (seriesUID.size() > 0)
		{
			std::cout << std::endl << "The directory: " << std::endl;
			std::cout << std::endl << argv[1] << std::endl << std::endl;
			std::cout << "Contains the following DICOM Series: ";
			std::cout << std::endl << std::endl;

			SeriesIdContainer::const_iterator seriesItr = seriesUID.begin();
			SeriesIdContainer::const_iterator seriesEnd = seriesUID.end();

			while (seriesItr != seriesEnd)
			{
				std::cout << seriesItr->c_str() << std::endl;
				++seriesItr;
			}

			std::string studyDescriptionOld;

			for (seriesItr = seriesUID.begin(); seriesItr != seriesEnd; seriesItr++)
			{
				std::string seriesIdentifier;
				seriesIdentifier = seriesItr->c_str();

				FileNamesContainer fileNames;
				fileNames = nameGenerator->GetFileNames(seriesIdentifier);

				std::cout << std::endl << std::endl;
				std::cout << "Now reading series: " << std::endl << std::endl;
				std::cout << seriesIdentifier << std::endl;
				std::cout << std::endl << std::endl;


				if (getModality(fileNames[0]) == "PT") // PET
				{
					readAndWrite<FloatImageType>(fileNames, seriesDirectory);
				}
				else // CT
				{
					readAndWrite<ShortImageType>(fileNames, seriesDirectory);
				}
			}
		}
	}
	catch (itk::ExceptionObject &ex)
	{
		std::cout << ex << std::endl;
		return EXIT_FAILURE;
	}
	

	return EXIT_SUCCESS;
}



<|start_filename|>Tools/FeatureExtraction/FeatureExtraction.cxx<|end_filename|>
#if defined(_MSC_VER)
# pragma warning ( disable : 4786 )
#endif

#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkResampleImageFilter.h>
#include <itkCastImageFilter.h>
#include <itkCommand.h>

#include <itkRegionOfInterestImageFilter.h>
#include <itkExtractImageFilter.h>

// Operations on Images

#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIterator.h>
#include <vnl/vnl_math.h>

// For Feature
#include <itkRescaleIntensityImageFilter.h>


#include <itkMinimumMaximumImageCalculator.h>
#include <itkStatisticsImageFilter.h>

#include <itkScalarImageToTextureFeaturesFilter.h>
// #include "itkScalarImageToRunLengthFeaturesFilter1.h"
#include <itkScalarImageToRunLengthFeaturesFilter.h>

#include <itkLabelGeometryImageFilter.h>
#include <itkConnectedComponentImageFilter.h>
#include <itkRelabelComponentImageFilter.h>
#include <itkLabelStatisticsImageFilter.h>

// For Label Object Representation
#include <itkLabelObject.h>
#include <itkLabelMap.h>
#include <itkShapeLabelObject.h>
#include <itkStatisticsLabelObject.h>
#include <itkBinaryImageToLabelMapFilter.h>
#include <itkLabelMapToLabelImageFilter.h>
#include <itkLabelImageToLabelMapFilter.h>
#include <itkShapeLabelMapFilter.h>
#include <itkStatisticsLabelMapFilter.h>
#include <itkBinaryImageToShapeLabelMapFilter.h>
#include <itkBinaryImageToStatisticsLabelMapFilter.h>

// For island removing
#include <itkBinaryShapeKeepNObjectsImageFilter.h>

#include <itkBinaryFillholeImageFilter.h>
#include <itkSliceBySliceImageFilter.h>
#include <itkSubtractImageFilter.h>
#include <itkAndImageFilter.h>


// For file operation
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>

// For threshold
#include <itkThresholdImageFilter.h>
#include <itkBinaryThresholdImageFilter.h>

#include "ITKUtils.h"

using namespace std;

unsigned int NumberOfBins = 64;

void
decodeCharactheristicGLCM(int x, char name[])
{
    switch (x) {
        case 0:
            strcpy(name, "Energy");
            break;
        case 1:
            strcpy(name, "Entropy");
            break;
        case 2:
            strcpy(name, "Correlation");
            break;
        case 3:
            strcpy(name, "InverseDifferenceMoment");
            break;
        case 4:
            strcpy(name, "Inertia");
            break;
        case 5:
            strcpy(name, "ClusterShade");
            break;
        case 6:
            strcpy(name, "ClusterProminence");
            break;
        case 7:
            strcpy(name, "HaralickCorrelation");
            break;
    }
}

void
decodeCharactheristicGLRM(int x, char name[])
{
    switch (x) {
        case 0:
            strcpy(name, "ShortRunEmphasis");
            break;
        case 1:
            strcpy(name, "LongRunEmphasis");
            break;
        case 2:
            strcpy(name, "GreyLevelNonuniformity");
            break;
        case 3:
            strcpy(name, "RunLengthNonuniformity");
            break;
        case 4:
            strcpy(name, "LowGreyLevelRunEmphasis");
            break;
        case 5:
            strcpy(name, "HighGreyLevelRunEmphasis");
            break;
        case 6:
            strcpy(name, "ShortRunLowGreyLevelEmphasis");
            break;
        case 7:
            strcpy(name, "ShortRunHighGreyLevelEmphasis");
            break;
        case 8:
            strcpy(name, "LongRunLowGreyLevelEmphasis");
            break;
        case 9:
            strcpy(name, "LongRunHighGreyLevelEmphasis");
            break;
    }
} // decodeCharactheristicGLRM

void
CalcIntensityFeatures(InputImageType::Pointer inputImage, MaskImageType::Pointer maskImage, InputPixelType &inputMin,
  InputPixelType &inputMax, ostream &outfile)
{
    using LabelStatisticsFilterType = itk::LabelStatisticsImageFilter<InputImageType, MaskImageType>;

    LabelStatisticsFilterType::Pointer LabelStatisticsfilter = LabelStatisticsFilterType::New();
    LabelStatisticsfilter->UseHistogramsOn();
    LabelStatisticsfilter->SetHistogramParameters(3000, -1024, 2000);

    LabelStatisticsfilter->SetInput(inputImage);
    LabelStatisticsfilter->SetLabelInput(maskImage);

    LabelStatisticsfilter->Update();
    cout << "inputImage features from the LabelStatisticsImageFilter:" << endl;
    // outfile << endl<< "inputImage features from the LabelStatisticsImageFilter:" << endl;

    cout << "Minimum = " << LabelStatisticsfilter->GetMinimum(maskValue) << endl;
    cout << "Maximum = " << LabelStatisticsfilter->GetMaximum(maskValue) << endl;
    cout << "Median  = " << LabelStatisticsfilter->GetMedian(maskValue) << endl;
    cout << "Mean    = " << LabelStatisticsfilter->GetMean(maskValue) << endl;
    cout << "Sigma   = " << LabelStatisticsfilter->GetSigma(maskValue) << endl;
    cout << "Variance = " << LabelStatisticsfilter->GetVariance(maskValue) << endl;
    cout << "Sum     = " << LabelStatisticsfilter->GetSum(maskValue) << endl;
    // cout << "Bounding box = " << LabelStatisticsfilter->GetBoundingBox( maskValue )  << endl;
    cout << "Region  = " << LabelStatisticsfilter->GetRegion(maskValue) << endl;
    // cout << "HistogramPointer  = " << LabelStatisticsfilter->GetHistogram( maskValue )   << endl;

    outfile << "Minimum = " << LabelStatisticsfilter->GetMinimum(maskValue) << endl;
    outfile << "Maximum = " << LabelStatisticsfilter->GetMaximum(maskValue) << endl;
    outfile << "Median  = " << LabelStatisticsfilter->GetMedian(maskValue) << endl;
    outfile << "Mean = " << LabelStatisticsfilter->GetMean(maskValue) << endl;
    outfile << "Sigma = " << LabelStatisticsfilter->GetSigma(maskValue) << endl;
    outfile << "Variance = " << LabelStatisticsfilter->GetVariance(maskValue) << endl;
    outfile << "Sum = " << LabelStatisticsfilter->GetSum(maskValue) << endl;
    // outfile << "Bounding box = " << LabelStatisticsfilter->GetBoundingBox( maskValue )  << endl;
    // outfile << "Region  = " << LabelStatisticsfilter->GetRegion( maskValue )   << endl;
    // cout << "HistogramPointer  = " << LabelStatisticsfilter->GetHistogram( maskValue )   << endl;


    inputMin = static_cast<InputPixelType>(LabelStatisticsfilter->GetMinimum(maskValue));
    inputMax = static_cast<InputPixelType>(LabelStatisticsfilter->GetMaximum(maskValue));
} // CalcIntensityFeatures

void
CalcGeometryFeatures(InputImageType::Pointer inputImage, MaskImageType::Pointer maskImage, ostream &outfile)
{
    SpacingType spacing = inputImage->GetSpacing();

    using LabelGeometryType = itk::LabelGeometryImageFilter<MaskImageType, InputImageType>;
    LabelGeometryType::Pointer labelGeometryFilter = LabelGeometryType::New();

    // Set up a connected components filter to label the binary objects.
    using ConnectedComponentType = itk::ConnectedComponentImageFilter<MaskImageType, MaskImageType>;
    ConnectedComponentType::Pointer connectedComponentFilter = ConnectedComponentType::New();
    connectedComponentFilter->SetInput(maskImage);

    // Relabel the components in order of size.
    using RelabelType = itk::RelabelComponentImageFilter<MaskImageType, MaskImageType>;
    RelabelType::Pointer relabeler = RelabelType::New();
    relabeler->SetInput(connectedComponentFilter->GetOutput());
    labelGeometryFilter->SetInput(relabeler->GetOutput());
    // labelGeometryFilter->SetInput(maskImage);

    // These generate optional outputs.
    labelGeometryFilter->CalculatePixelIndicesOn();
    labelGeometryFilter->CalculateOrientedBoundingBoxOn();
    labelGeometryFilter->CalculateOrientedLabelRegionsOn();
    labelGeometryFilter->CalculateOrientedIntensityRegionsOn();

    labelGeometryFilter->ReleaseDataFlagOn();

    labelGeometryFilter->SetIntensityInput(inputImage); //


    labelGeometryFilter->Update();
    cout << "Image features from the LabelGeometryImageFilter:" << endl;
    // outfile << "inputImage features from the LabelGeometryImageFilter:" << endl;

    outfile << "Volume = " << labelGeometryFilter->GetVolume(maskValue) * spacing[0] * spacing[1] * spacing[2] / 1000
            << endl;
    outfile << "IntegratedIntensity = " << labelGeometryFilter->GetIntegratedIntensity(maskValue) << endl;
    outfile << "Centroid = " << labelGeometryFilter->GetCentroid(maskValue) << endl;
    outfile << "WeightedCentroid = " << labelGeometryFilter->GetWeightedCentroid(maskValue) << endl;
    outfile << "AxesLength = " << labelGeometryFilter->GetAxesLength(maskValue) << endl;
    outfile << "MajorAxisLength = " << labelGeometryFilter->GetMajorAxisLength(maskValue) << endl;
    outfile << "MinorAxisLength = " << labelGeometryFilter->GetMinorAxisLength(maskValue) << endl;
    outfile << "Eccentricity = " << labelGeometryFilter->GetEccentricity(maskValue) << endl;
    outfile << "Elongation = " << labelGeometryFilter->GetElongation(maskValue) << endl;
    outfile << "Orientation = " << labelGeometryFilter->GetOrientation(maskValue) << endl;
    // outfile << "Boundingbox = " << labelGeometryFilter->GetBoundingBox(maskValue) << endl;

    outfile << "BoundingBoxVolume = " << labelGeometryFilter->GetBoundingBoxVolume(maskValue) * spacing[0]
        * spacing[1] * spacing[2] / 1000 << endl;
    outfile << "BoundingBoxSize = " << labelGeometryFilter->GetBoundingBoxSize(maskValue) << endl;
    outfile << "OrientedBoundingBoxVolume = " << labelGeometryFilter->GetOrientedBoundingBoxVolume(maskValue)
        * spacing[0] * spacing[1] * spacing[2] / 1000 << endl;
    outfile << "OrientedBoundingBoxSize = " << labelGeometryFilter->GetOrientedBoundingBoxSize(maskValue) << endl;

    // outfile << "Eigenvalues = " << labelGeometryFilter->GetEigenvalues(maskValue) << endl;
    outfile << "Eigenvalues= [" << labelGeometryFilter->GetEigenvalues(maskValue)[0] << ", "
            << labelGeometryFilter->GetEigenvalues(maskValue)[1] << ", "
            << labelGeometryFilter->GetEigenvalues(maskValue)[2] << "]" << endl;
    outfile << "Eigenvectors = " << labelGeometryFilter->GetEigenvectors(maskValue) << endl;
    outfile << "RotationMatrix = " << labelGeometryFilter->GetRotationMatrix(maskValue) << endl;
    // outfile << "\n\n";
} // CalcGeometryFeatures

void
CalcShapeAndIntensityFeatures(InputImageType::Pointer inputImage, MaskImageType::Pointer maskImage, ostream &outfile)
{
    using LabelObjectTypeS = itk::StatisticsLabelObject<LabelPixelType, Dimension>;
    using LabelMapTypeS = itk::LabelMap<LabelObjectTypeS>;

    // converting binary image to Statistics label map
    using I2LSType = itk::BinaryImageToStatisticsLabelMapFilter<MaskImageType, InputImageType, LabelMapTypeS>;
    I2LSType::Pointer i2ls = I2LSType::New();
    // itk::SimpleFilterWatcher watcher1( i2ls );

    i2ls->FullyConnectedOn();
    int inputForegroundValue  = 1;
    bool computeFeretDiameter = 1;
    unsigned int outputBackgroundValue = 0;
    bool computePerimeter     = 1;
    bool computeHistogram     = 1;
    unsigned int numberOfBins = NumberOfBins;
    i2ls->SetInputForegroundValue(inputForegroundValue);
    i2ls->SetOutputBackgroundValue(outputBackgroundValue);
    i2ls->SetComputeFeretDiameter(computeFeretDiameter);
    i2ls->ComputeFeretDiameterOff();
    i2ls->ComputeFeretDiameterOn();
    i2ls->SetComputePerimeter(computePerimeter);
    i2ls->ComputePerimeterOn();
    i2ls->SetComputeHistogram(computeHistogram);
    i2ls->ComputeHistogramOn();
    i2ls->SetNumberOfBins(numberOfBins);


    LabelMapTypeS::Pointer labelMapS;
    LabelObjectTypeS * labelObjectS;
    i2ls->SetFeatureImage(inputImage);
    i2ls->SetInput(maskImage);

    i2ls->Update();
    cout << "inputImage features from the BinaryImageToStatisticsLabelMapFilter:" << endl;
    // outfile << "inputImage features from the BinaryImageToStatisticsLabelMapFilter:" << endl;


    labelMapS    = i2ls->GetOutput();
    labelObjectS = labelMapS->GetLabelObject(maskValue);


    //    cout << "BoundingBox = " << labelObjectS->GetBoundingBox() << endl;
    cout << "CenterOfGravity = " << labelObjectS->GetCenterOfGravity() << endl;
    cout << "Centroid = " << labelObjectS->GetCentroid() << endl;
    cout << "Elongation = " << labelObjectS->GetElongation() << endl;
    //  cout << "EquivalentEllipsoidDiameter = " << labelObjectS->GetEquivalentEllipsoidDiameter() << endl;
    //  cout << "EquivalentSphericalPerimeter = " << labelObjectS->GetEquivalentSphericalPerimeter() << endl;
    //       cout << "EquivalentSphericalRadius = " << labelObjectS->GetEquivalentSphericalRadius() << endl;
    cout << "FeretDiameter = " << labelObjectS->GetFeretDiameter() << endl;
    // cout << "Histogram = " << labelObjectS->GetHistogram() << endl;
    //  cout << "Index = " << labelObjectS->GetIndex() << endl;
    cout << "Kurtosis = " << labelObjectS->GetKurtosis() << endl;
    cout << "Flatness = " << labelObjectS->GetFlatness() << endl;
    cout << "Label = " << labelObjectS->GetLabel() << endl;
    //  cout << "Line = " << labelObjectS->GetLine() << endl;
    //        cout << "LineContainer = " << labelObjectS->GetLineContainer() << endl;
    cout << "Maximum = " << labelObjectS->GetMaximum() << endl;
    // cout << "MaximumIndex = " << labelObjectS->GetMaximumIndex() << endl;
    cout << "Mean = " << labelObjectS->GetMean() << endl;
    cout << "Median = " << labelObjectS->GetMedian() << endl;
    cout << "Minimum = " << labelObjectS->GetMinimum() << endl;
    // cout << "MinimumIndex = " << labelObjectS->GetMinimumIndex() << endl;
    cout << "NameOfClass = " << labelObjectS->GetNameOfClass() << endl;
    cout << "NumberOfLines = " << labelObjectS->GetNumberOfLines() << endl;
    //        cout << "NumberOfPixels = " << labelObjectS->GetNumberOfPixels() << endl;
    //  cout << "NumberOfPixelsOnBorder = " << labelObjectS->GetNumberOfPixelsOnBorder() << endl;
    cout << "Perimeter = " << labelObjectS->GetPerimeter() << endl;
    //  cout << "PerimeterOnBorder = " << labelObjectS->GetPerimeterOnBorder() << endl;
    //        cout << "PerimeterOnBorderRatio = " << labelObjectS->GetPerimeterOnBorderRatio() << endl;
    cout << "PhysicalSize = " << labelObjectS->GetPhysicalSize() << endl;
    cout << "PrincipalAxes = " << labelObjectS->GetPrincipalAxes() << endl;
    cout << "PrincipalMoments = " << labelObjectS->GetPrincipalMoments() << endl;
    cout << "ReferenceCount = " << labelObjectS->GetReferenceCount() << endl;
    cout << "RegionElongation = " << labelObjectS->GetElongation() << endl; // RegionElongation is disappeared
    cout << "Roundness = " << labelObjectS->GetRoundness() << endl;
    cout << "SizeRegionRatio = " << 0 /*labelObjectS->GetSizeRegionRatio()*/ << endl; // This property is no longer provide
    cout << "Skewness = " << labelObjectS->GetSkewness() << endl;
    cout << "StandardDeviation = " << labelObjectS->GetStandardDeviation() << endl;
    cout << "Sum = " << labelObjectS->GetSum() << endl;
    cout << "Variance = " << labelObjectS->GetVariance() << endl;
    //  cout << "WeightedElongation = " << labelObjectS->GetWeightedElongation() << endl;
    //        cout << "WeightedFlatness = " << labelObjectS->GetWeightedFlatness() << endl;
    //  cout << "WeightedPrincipalAxes = " << labelObjectS->GetWeightedPrincipalAxes() << endl;
    //        cout << "WeightedPrincipalMoments = " << labelObjectS->GetWeightedPrincipalMoments() << endl;

    // outfile << "BoundingBox = " << labelObjectS->GetBoundingBox() << endl;
    outfile << "CenterOfGravity = " << labelObjectS->GetCenterOfGravity() << endl;
    outfile << "Centroid = " << labelObjectS->GetCentroid() << endl;
    outfile << "Elongation = " << labelObjectS->GetElongation() << endl;
    outfile << "EquivalentEllipsoidDiameter = " << labelObjectS->GetEquivalentEllipsoidDiameter() << endl;
    outfile << "EquivalentSphericalPerimeter = " << labelObjectS->GetEquivalentSphericalPerimeter() << endl;
    outfile << "EquivalentSphericalRadius = " << labelObjectS->GetEquivalentSphericalRadius() << endl;
    outfile << "FeretDiameter = " << labelObjectS->GetFeretDiameter() << endl;
    // outfile << "Histogram = " << labelObjectS->GetHistogram() << endl;
    // outfile << "Index = " << labelObjectS->GetIndex() << endl;
    outfile << "Kurtosis = " << labelObjectS->GetKurtosis() << endl;
    outfile << "Flatness = " << labelObjectS->GetFlatness() << endl;
    // outfile << "Label = " << labelObjectS->GetLabel() << endl;
    // outfile << "Line = " << labelObjectS->GetLine() << endl;
    // outfile << "LineContainer = " << labelObjectS->GetLineContainer() << endl;
    outfile << "Maximum = " << labelObjectS->GetMaximum() << endl;
    // outfile << "Maximum Index = " << labelObjectS->GetMaximumIndex() << endl;
    outfile << "Mean = " << labelObjectS->GetMean() << endl;
    outfile << "Median = " << labelObjectS->GetMedian() << endl;
    outfile << "Minimum = " << labelObjectS->GetMinimum() << endl;
    // outfile << "MinimumIndex = " << labelObjectS->GetMinimumIndex() << endl;
    // outfile << "NameOfClass = " << labelObjectS->GetNameOfClass() << endl;
    outfile << "NumberOfLines = " << labelObjectS->GetNumberOfLines() << endl;
    outfile << "NumberOfPixels = " << labelObjectS->GetNumberOfPixels() << endl;
    outfile << "NumberOfPixelsOnBorder = " << labelObjectS->GetNumberOfPixelsOnBorder() << endl;
    outfile << "Perimeter = " << labelObjectS->GetPerimeter() << endl;
    outfile << "PerimeterOnBorder = " << labelObjectS->GetPerimeterOnBorder() << endl;
    outfile << "PerimeterOnBorderRatio = " << labelObjectS->GetPerimeterOnBorderRatio() << endl;
    outfile << "PhysicalSize = " << labelObjectS->GetPhysicalSize() << endl;
    outfile << "PrincipalAxes = " << labelObjectS->GetPrincipalAxes() << endl;
    outfile << "PrincipalMoments = " << labelObjectS->GetPrincipalMoments() << endl;
    outfile << "ReferenceCount = " << labelObjectS->GetReferenceCount() << endl;
    outfile << "RegionElongation = " << labelObjectS->GetElongation() << endl; // deleted
    outfile << "Roundness = " << labelObjectS->GetRoundness() << endl;
    // outfile << "SizeRegionRatio = " << labelObjectS->GetSizeRegionRatio() << endl; // deleted
    outfile << "Skewness = " << labelObjectS->GetSkewness() << endl;
    outfile << "StandardDeviation = " << labelObjectS->GetStandardDeviation() << endl;
    outfile << "Sum = " << labelObjectS->GetSum() << endl;
    outfile << "Variance = " << labelObjectS->GetVariance() << endl;
    outfile << "WeightedElongation = " << labelObjectS->GetWeightedElongation() << endl;
    outfile << "WeightedFlatness = " << labelObjectS->GetWeightedFlatness() << endl;
    outfile << "WeightedPrincipalAxes = " << labelObjectS->GetWeightedPrincipalAxes() << endl;
    outfile << "WeightedPrincipalMoments = " << labelObjectS->GetWeightedPrincipalMoments() << endl;
} // CalcShapeAndIntensityFeatures

void
CalcShapeAndIntensityFeatures(InputImage2DType::Pointer inputImage, MaskImage2DType::Pointer maskImage,
  ostream &outfile)
{
    using LabelObject2DTypeS = itk::StatisticsLabelObject<LabelPixelType, Dimension - 1>;
    using LabelMap2DTypeS = itk::LabelMap<LabelObject2DTypeS>;

    // converting binary image to Statistics label map
    using I2LS2DType = itk::BinaryImageToStatisticsLabelMapFilter<MaskImage2DType, InputImage2DType, LabelMap2DTypeS>;
    I2LS2DType::Pointer i2ls = I2LS2DType::New();
    // itk::SimpleFilterWatcher watcher1( i2ls );

    i2ls->FullyConnectedOn();
    int inputForegroundValue  = 1;
    bool computeFeretDiameter = 1;
    unsigned int outputBackgroundValue = 0;
    bool computePerimeter     = 1;
    bool computeHistogram     = 1;
    unsigned int numberOfBins = NumberOfBins;
    i2ls->SetInputForegroundValue(inputForegroundValue);
    i2ls->SetOutputBackgroundValue(outputBackgroundValue);
    i2ls->SetComputeFeretDiameter(computeFeretDiameter);
    i2ls->ComputeFeretDiameterOff();
    i2ls->ComputeFeretDiameterOn();
    i2ls->SetComputePerimeter(computePerimeter);
    i2ls->ComputePerimeterOn();
    i2ls->SetComputeHistogram(computeHistogram);
    i2ls->ComputeHistogramOn();
    i2ls->SetNumberOfBins(numberOfBins);


    LabelMap2DTypeS::Pointer labelMapS;
    LabelObject2DTypeS * labelObjectS;
    i2ls->SetFeatureImage(inputImage);
    i2ls->SetInput(maskImage);

    i2ls->Update();
    cout << "inputImage features from the BinaryImageToStatisticsLabelMapFilter:" << endl;
    // outfile << "inputImage features from the BinaryImageToStatisticsLabelMapFilter:" << endl;


    labelMapS    = i2ls->GetOutput();
    labelObjectS = labelMapS->GetLabelObject(maskValue);


    //    cout << "2D BoundingBox = " << labelObjectS->GetBoundingBox() << endl;
    cout << "2D CenterOfGravity = " << labelObjectS->GetCenterOfGravity() << endl;
    cout << "2D Centroid = " << labelObjectS->GetCentroid() << endl;
    cout << "2D Elongation = " << labelObjectS->GetElongation() << endl;
    cout << "2D EquivalentEllipsoidDiameter = " << labelObjectS->GetEquivalentEllipsoidDiameter() << endl;
    cout << "2D EquivalentSphericalPerimeter = " << labelObjectS->GetEquivalentSphericalPerimeter() << endl;
    cout << "2D EquivalentSphericalRadius = " << labelObjectS->GetEquivalentSphericalRadius() << endl;
    cout << "2D FeretDiameter = " << labelObjectS->GetFeretDiameter() << endl;
    // cout << "2D Histogram = " << labelObjectS->GetHistogram() << endl;
    //  cout << "2D Index = " << labelObjectS->GetIndex() << endl;
    cout << "2D Kurtosis = " << labelObjectS->GetKurtosis() << endl;
    cout << "2D Flatness = " << labelObjectS->GetFlatness() << endl;
    // cout << "2D Label = " << labelObjectS->GetLabel() << endl;
    //  cout << "2D Line = " << labelObjectS->GetLine() << endl;
    //        cout << "2D LineContainer = " << labelObjectS->GetLineContainer() << endl;
    cout << "2D Maximum = " << labelObjectS->GetMaximum() << endl;
    // cout << "2D MaximumIndex = " << labelObjectS->GetMaximumIndex() << endl;
    cout << "2D Mean = " << labelObjectS->GetMean() << endl;
    cout << "2D Median = " << labelObjectS->GetMedian() << endl;
    cout << "2D Minimum = " << labelObjectS->GetMinimum() << endl;
    // cout << "2D MinimumIndex = " << labelObjectS->GetMinimumIndex() << endl;
    cout << "2D NameOfClass = " << labelObjectS->GetNameOfClass() << endl;
    cout << "2D NumberOfLines = " << labelObjectS->GetNumberOfLines() << endl;
    //        cout << "2D NumberOfPixels = " << labelObjectS->GetNumberOfPixels() << endl;
    //  cout << "2D NumberOfPixelsOnBorder = " << labelObjectS->GetNumberOfPixelsOnBorder() << endl;
    cout << "2D Perimeter = " << labelObjectS->GetPerimeter() << endl;
    cout << "2D PerimeterOnBorder = " << labelObjectS->GetPerimeterOnBorder() << endl;
    cout << "2D PerimeterOnBorderRatio = " << labelObjectS->GetPerimeterOnBorderRatio() << endl;
    cout << "2D PhysicalSize = " << labelObjectS->GetPhysicalSize() << endl;
    cout << "2D PrincipalAxes = " << labelObjectS->GetPrincipalAxes() << endl;
    cout << "2D PrincipalMoments = " << labelObjectS->GetPrincipalMoments() << endl;
    cout << "2D ReferenceCount = " << labelObjectS->GetReferenceCount() << endl;
    cout << "2D RegionElongation = " << labelObjectS->GetElongation() << endl; // RegionElongation is disappeared
    cout << "2D Roundness = " << labelObjectS->GetRoundness() << endl;
    cout << "2D SizeRegionRatio = " << 0 /*labelObjectS->GetSizeRegionRatio()*/ << endl; // This property is no longer provide
    cout << "2D Skewness = " << labelObjectS->GetSkewness() << endl;
    cout << "2D StandardDeviation = " << labelObjectS->GetStandardDeviation() << endl;
    cout << "2D Sum = " << labelObjectS->GetSum() << endl;
    cout << "2D Variance = " << labelObjectS->GetVariance() << endl;
    //  cout << "2D WeightedElongation = " << labelObjectS->GetWeightedElongation() << endl;
    //        cout << "2D WeightedFlatness = " << labelObjectS->GetWeightedFlatness() << endl;
    //  cout << "2D WeightedPrincipalAxes = " << labelObjectS->GetWeightedPrincipalAxes() << endl;
    //        cout << "2D WeightedPrincipalMoments = " << labelObjectS->GetWeightedPrincipalMoments() << endl;

    // outfile << "2DBoundingBox = " << labelObjectS->GetBoundingBox() << endl;
    outfile << "2DCenterOfGravity = " << labelObjectS->GetCenterOfGravity() << endl;
    outfile << "2DCentroid = " << labelObjectS->GetCentroid() << endl;
    outfile << "2DElongation = " << labelObjectS->GetElongation() << endl;
    outfile << "2DEquivalentEllipsoidDiameter = " << labelObjectS->GetEquivalentEllipsoidDiameter() << endl;
    outfile << "2DEquivalentSphericalPerimeter = " << labelObjectS->GetEquivalentSphericalPerimeter() << endl;
    outfile << "2DEquivalentSphericalRadius = " << labelObjectS->GetEquivalentSphericalRadius() << endl;
    outfile << "2DFeretDiameter = " << labelObjectS->GetFeretDiameter() << endl;
    // outfile << "2DHistogram = " << labelObjectS->GetHistogram() << endl;
    // outfile << "2DIndex = " << labelObjectS->GetIndex() << endl;
    outfile << "2DKurtosis = " << labelObjectS->GetKurtosis() << endl;
    outfile << "2DFlatness = " << labelObjectS->GetFlatness() << endl;
    // outfile << "2DLabel = " << labelObjectS->GetLabel() << endl;
    // outfile << "2DLine = " << labelObjectS->GetLine() << endl;
    // outfile << "2DLineContainer = " << labelObjectS->GetLineContainer() << endl;
    outfile << "2DMaximum = " << labelObjectS->GetMaximum() << endl;
    // outfile << "2DMaximum Index = " << labelObjectS->GetMaximumIndex() << endl;
    outfile << "2DMean = " << labelObjectS->GetMean() << endl;
    outfile << "2DMedian = " << labelObjectS->GetMedian() << endl;
    outfile << "2DMinimum = " << labelObjectS->GetMinimum() << endl;
    // outfile << "2DMinimumIndex = " << labelObjectS->GetMinimumIndex() << endl;
    // outfile << "2DNameOfClass = " << labelObjectS->GetNameOfClass() << endl;
    outfile << "2DNumberOfLines = " << labelObjectS->GetNumberOfLines() << endl;
    outfile << "2DNumberOfPixels = " << labelObjectS->GetNumberOfPixels() << endl;
    outfile << "2DNumberOfPixelsOnBorder = " << labelObjectS->GetNumberOfPixelsOnBorder() << endl;
    outfile << "2DPerimeter = " << labelObjectS->GetPerimeter() << endl;
    outfile << "2DPerimeterOnBorder = " << labelObjectS->GetPerimeterOnBorder() << endl;
    outfile << "2DPerimeterOnBorderRatio = " << labelObjectS->GetPerimeterOnBorderRatio() << endl;
    outfile << "2DPhysicalSize = " << labelObjectS->GetPhysicalSize() << endl;
    outfile << "2DPrincipalAxes = " << labelObjectS->GetPrincipalAxes() << endl;
    outfile << "2DPrincipalMoments = " << labelObjectS->GetPrincipalMoments() << endl;
    outfile << "2DReferenceCount = " << labelObjectS->GetReferenceCount() << endl;
    outfile << "2DRegionElongation = " << labelObjectS->GetElongation() << endl; // deleted
    outfile << "2DRoundness = " << labelObjectS->GetRoundness() << endl;
    // outfile << "2DSizeRegionRatio = " << labelObjectS->GetSizeRegionRatio() << endl; // deleted
    outfile << "2DSkewness = " << labelObjectS->GetSkewness() << endl;
    outfile << "2DStandardDeviation = " << labelObjectS->GetStandardDeviation() << endl;
    outfile << "2DSum = " << labelObjectS->GetSum() << endl;
    outfile << "2DVariance = " << labelObjectS->GetVariance() << endl;
    outfile << "2DWeightedElongation = " << labelObjectS->GetWeightedElongation() << endl;
    outfile << "2DWeightedFlatness = " << labelObjectS->GetWeightedFlatness() << endl;
    outfile << "2DWeightedPrincipalAxes = " << labelObjectS->GetWeightedPrincipalAxes() << endl;
    outfile << "2DWeightedPrincipalMoments = " << labelObjectS->GetWeightedPrincipalMoments() << endl;
} // CalcShapeAndIntensityFeatures

void
CalcGlcmFeatures(InputImageType::Pointer inputImage, InputImageType::Pointer maskImage, InputPixelType inputMin,
  InputPixelType inputMax, ostream &outfile)
{
    using TextureFilterType = itk::Statistics::ScalarImageToTextureFeaturesFilter<InputImageType>;
    TextureFilterType::Pointer glcm = TextureFilterType::New();

    using TextureFeaturesFilterType = TextureFilterType::TextureFeaturesFilterType;
    using TextureFeatureEnums = itk::Statistics::HistogramToTextureFeaturesFilterEnums::TextureFeature;
    TextureFilterType::FeatureNameVectorPointer requestedFeatures =
      TextureFilterType::FeatureNameVector::New();

    requestedFeatures->push_back(static_cast<uint8_t>(TextureFeatureEnums::Energy));
    requestedFeatures->push_back(static_cast<uint8_t>(TextureFeatureEnums::Entropy));
    requestedFeatures->push_back(static_cast<uint8_t>(TextureFeatureEnums::Correlation));
    requestedFeatures->push_back(static_cast<uint8_t>(TextureFeatureEnums::InverseDifferenceMoment));
    requestedFeatures->push_back(static_cast<uint8_t>(TextureFeatureEnums::Inertia));
    requestedFeatures->push_back(static_cast<uint8_t>(TextureFeatureEnums::ClusterShade));
    requestedFeatures->push_back(static_cast<uint8_t>(TextureFeatureEnums::ClusterProminence));
    requestedFeatures->push_back(static_cast<uint8_t>(TextureFeatureEnums::HaralickCorrelation));

    glcm->SetRequestedFeatures(requestedFeatures);
    // glcm->SetNormalizeOn();
    // glcm->Print( cout );


    glcm->SetInput(inputImage);

    glcm->SetMaskImage(maskImage);
    // glcm->FastCalculationsOn();

    glcm->SetNumberOfBinsPerAxis(NumberOfBins);
    glcm->SetPixelValueMinMax(inputMin, inputMax);

    glcm->FastCalculationsOff();
    glcm->Update();
    cout << "The inputImage features from the ScalarImageToTextureFeaturesFilter:" << endl;
    // outfile<< endl<< endl << "The inputImage features from the ScalarImageToTextureFeaturesFilter:" << endl;


    using GLCMFilterType = itk::Statistics::ScalarImageToCooccurrenceMatrixFilter<InputImageType>;
    using HistogramType = GLCMFilterType::HistogramType;
    using Hist2FeaturesType = itk::Statistics::HistogramToTextureFeaturesFilter<HistogramType>;

    TextureFilterType::OffsetVector::ConstPointer offsets = glcm->GetOffsets();
    for (TextureFilterType::OffsetVector::ConstIterator offset = offsets->Begin(); offset != offsets->End(); offset++) {
        GLCMFilterType::Pointer glcmGenerator = GLCMFilterType::New();
        glcmGenerator->SetOffset(offset->Value());
        glcmGenerator->SetNumberOfBinsPerAxis(NumberOfBins);    // reasonable number of bins
        glcmGenerator->SetPixelValueMinMax(inputMin, inputMax); // for input UCHAR pixel type
        Hist2FeaturesType::Pointer featureCalc = Hist2FeaturesType::New();

        glcmGenerator->SetInput(inputImage);
        glcmGenerator->SetMaskImage(maskImage);
        glcmGenerator->Update();
        featureCalc->SetInput(glcmGenerator->GetOutput());
        featureCalc->Update();

        cout << offset->Value() << "\t" << featureCalc->GetFeature(Hist2FeaturesType::InverseDifferenceMoment) << endl;
    }


    int fsize = glcm->GetFeatureMeans()->Size();
    for (int x = 0; x < fsize; x++) {
        char charactheristic[100];
        decodeCharactheristicGLCM(glcm->GetRequestedFeatures()->GetElement(x), charactheristic);

        cout << "MeanOf" << charactheristic << "= " << glcm->GetFeatureMeans()->GetElement(x) << endl;

        outfile << "MeanOf" << charactheristic << "= " << glcm->GetFeatureMeans()->GetElement(x) << endl;
    }
    for (int x = 0; x < fsize; x++) {
        char charactheristic[100];
        decodeCharactheristicGLCM(glcm->GetRequestedFeatures()->GetElement(x), charactheristic);

        cout << "StandardDeviationOf" << charactheristic << "= "
             << glcm->GetFeatureStandardDeviations()->GetElement(x) << endl;

        outfile << "StandardDeviationOf" << charactheristic << "= " << glcm->GetFeatureStandardDeviations()->GetElement(
            x) << endl;
    }
} // CalcGlcmFeatures

void
CalcGlrmFeatures(InputImageType::Pointer inputImage, InputImageType::Pointer maskImage, InputPixelType inputMin,
  InputPixelType inputMax, double distMax, ostream &outfile)
{
    using RunLengthFilterType = itk::Statistics::ScalarImageToRunLengthFeaturesFilter<InputImageType>;
    RunLengthFilterType::Pointer glrm = RunLengthFilterType::New();

    using RunLengthFeaturesFilterType = RunLengthFilterType::RunLengthFeaturesFilterType;
    using RunLengthFeatureEnums = itk::Statistics::HistogramToRunLengthFeaturesFilterEnums::RunLengthFeature;
    RunLengthFilterType::FeatureNameVectorPointer requestedFeatures =
      RunLengthFilterType::FeatureNameVector::New();

    requestedFeatures->push_back(static_cast<uint8_t>(RunLengthFeatureEnums::ShortRunEmphasis));
    requestedFeatures->push_back(static_cast<uint8_t>(RunLengthFeatureEnums::LongRunEmphasis));
    requestedFeatures->push_back(static_cast<uint8_t>(RunLengthFeatureEnums::GreyLevelNonuniformity));
    requestedFeatures->push_back(static_cast<uint8_t>(RunLengthFeatureEnums::RunLengthNonuniformity));
    requestedFeatures->push_back(static_cast<uint8_t>(RunLengthFeatureEnums::LowGreyLevelRunEmphasis));
    requestedFeatures->push_back(static_cast<uint8_t>(RunLengthFeatureEnums::HighGreyLevelRunEmphasis));
    requestedFeatures->push_back(static_cast<uint8_t>(RunLengthFeatureEnums::ShortRunLowGreyLevelEmphasis));
    requestedFeatures->push_back(static_cast<uint8_t>(RunLengthFeatureEnums::ShortRunHighGreyLevelEmphasis));
    requestedFeatures->push_back(static_cast<uint8_t>(RunLengthFeatureEnums::LongRunLowGreyLevelEmphasis));
    requestedFeatures->push_back(static_cast<uint8_t>(RunLengthFeatureEnums::LongRunHighGreyLevelEmphasis));

    glrm->SetRequestedFeatures(requestedFeatures);


    glrm->SetInput(inputImage);
    // glrm->SetInput(maskImage);
    glrm->SetMaskImage(maskImage);
    // glrm->FastCalculationsOn();

    glrm->SetNumberOfBinsPerAxis(NumberOfBins);
    glrm->SetPixelValueMinMax(inputMin, (inputMax - inputMin) * 2 + inputMin);
    glrm->SetDistanceValueMinMax(0, distMax * 2 / 3);

    glrm->FastCalculationsOff();
    glrm->Update();
    cout << "The inputImage features from the ScalarImageToRunLengthFeaturesFilter:" << endl;
    // outfile<< endl<< endl << "The inputImage features from the ScalarImageToTextureFeaturesFilter:" << endl;


    int size = glrm->GetFeatureMeans()->Size();
    for (int x = 0; x < size; x++) {
        char charactheristic[100];
        decodeCharactheristicGLRM(glrm->GetRequestedFeatures()->GetElement(x), charactheristic);

        cout << "MeanOf" << charactheristic << "= " << glrm->GetFeatureMeans()->GetElement(x) << endl;

        outfile << "MeanOf" << charactheristic << "= " << glrm->GetFeatureMeans()->GetElement(x) << endl;
    }
    for (int x = 0; x < size; x++) {
        char charactheristic[100];
        decodeCharactheristicGLRM(glrm->GetRequestedFeatures()->GetElement(x), charactheristic);

        cout << "StandardDeviationOf" << charactheristic << "= "
             << glrm->GetFeatureStandardDeviations()->GetElement(x) << endl;

        outfile << "StandardDeviationOf" << charactheristic << "= " << glrm->GetFeatureStandardDeviations()->GetElement(
            x) << endl;
    }
} // CalcGlrmFeatures

MaskImageType::Pointer
StentRemoval(InputImageType::Pointer inputImage, MaskImageType::Pointer maskImage, InputPixelType stentThreshold = 500)
{
    StructuringElementType structuringElement1;

    structuringElement1.SetRadius(1);
    structuringElement1.CreateStructuringElement();

    StructuringElementType structuringElement3;
    structuringElement3.SetRadius(3);
    structuringElement3.CreateStructuringElement();

    StructuringElementType structuringElement5;
    structuringElement5.SetRadius(5);
    structuringElement5.CreateStructuringElement();

    StructuringElementType structuringElement10;
    structuringElement10.SetRadius(10);
    structuringElement10.CreateStructuringElement();

    #if 0
    ErodeFilterType::Pointer erodeFilter = ErodeFilterType::New();
    erodeFilter->SetInput(maskImage);
    erodeFilter->SetKernel(structuringElement3);
    erodeFilter->SetErodeValue(maskValue);
    erodeFilter->Update();

    maskImage = erodeFilter->GetOutput();
    #endif


    using BinaryThresholdfilterType = itk::BinaryThresholdImageFilter<InputImageType, MaskImageType>;

    BinaryThresholdfilterType::Pointer binaryThresholdfilter = BinaryThresholdfilterType::New();

    binaryThresholdfilter->SetInput(inputImage);
    binaryThresholdfilter->SetInsideValue(1);
    binaryThresholdfilter->SetOutsideValue(0);
    binaryThresholdfilter->SetLowerThreshold(500); // <- Main paraneter

    #if 0
    using LabelOpeningType = itk::BinaryShapeKeepNObjectsImageFilter<MaskImageType>;

    LabelOpeningType::Pointer opening = LabelOpeningType::New();
    opening->SetInput(binaryThresholdfilter->GetOutput());
    opening->SetBackgroundValue(0);
    opening->SetForegroundValue(maskValue);
    opening->SetNumberOfObjects(1);
    // opening->SetReverseOrdering( false );
    // opening->SetAttribute(100); //  extract the largest object  NUMBER_OF_PIXELS
    opening->Update();
    #endif // if 0

    DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
    dilateFilter->SetInput(binaryThresholdfilter->GetOutput());
    dilateFilter->SetKernel(structuringElement3);
    dilateFilter->SetDilateValue(maskValue);
    dilateFilter->Update();

    using MaskImage2DType = itk::Image<MaskPixelType, Dimension - 1>;

    using SliceBySliceFilterType = itk::SliceBySliceImageFilter<MaskImageType, MaskImageType>;
    SliceBySliceFilterType::Pointer sliceBySliceFilter = SliceBySliceFilterType::New();

    using I2LType = itk::BinaryFillholeImageFilter<MaskImage2DType>;
    I2LType::Pointer reconstruction = I2LType::New();
    reconstruction->SetFullyConnected(true);
    reconstruction->SetForegroundValue(maskValue);

    sliceBySliceFilter->SetInput(dilateFilter->GetOutput());
    sliceBySliceFilter->SetFilter(reconstruction);

    sliceBySliceFilter->Update();


    // remove stent area
    using SubtractImageFilterType = itk::SubtractImageFilter<MaskImageType, MaskImageType>;
    SubtractImageFilterType::Pointer subtractFilter = SubtractImageFilterType::New();
    subtractFilter->SetInput1(maskImage);
    subtractFilter->SetInput2(sliceBySliceFilter->GetOutput());

    using AndImageFilterType = itk::AndImageFilter<MaskImageType>;
    AndImageFilterType::Pointer andFilter = AndImageFilterType::New();
    andFilter->SetInput1(maskImage);
    andFilter->SetInput2(subtractFilter->GetOutput());
    andFilter->Update();
    maskImage = andFilter->GetOutput();

    return maskImage;
} // StentRemoval

int
ExtractLargestAreaSlice(InputImageType::Pointer inputImage, MaskImageType::Pointer maskImage,
  InputImage2DType::Pointer& input2DImage, MaskImage2DType::Pointer& mask2DImage)
{
    MaskImageType::RegionType maskRegion = maskImage->GetLargestPossibleRegion();
    MaskImageType::IndexType start       = maskRegion.GetIndex();
    MaskImageType::SizeType size         = maskRegion.GetSize();

    const unsigned int numberofSlices = size[2];
    size[2] = 0;


    using FilterType = itk::ExtractImageFilter<MaskImageType, MaskImage2DType>;
    FilterType::Pointer filter = FilterType::New();
    filter->SetDirectionCollapseToSubmatrix();

    float maxArea         = 0;
    unsigned int maxSlice = 0;
    for (unsigned int i = 0; i < numberofSlices; i++) {
        start[2] = i;

        MaskImageType::RegionType desiredRegion;
        desiredRegion.SetSize(size);
        desiredRegion.SetIndex(start);

        filter->SetExtractionRegion(desiredRegion);
        filter->SetInput(maskImage);
        filter->Update();

        using maskImage2DIteratorType = itk::ImageRegionConstIterator<MaskImage2DType>;
        maskImage2DIteratorType maskImage2Dit(filter->GetOutput(), filter->GetOutput()->GetRequestedRegion());

        float curArea = 0;
        maskImage2Dit.GoToBegin();
        while (!maskImage2Dit.IsAtEnd()) {
            if (maskImage2Dit.Get() == maskValue) ++curArea;
            ++maskImage2Dit;
        }
        if (curArea > maxArea) {
            maxArea  = curArea;
            maxSlice = i;
        }
        cout << i << ": max area " << maxArea << ", current area " << curArea << endl;
    }


    // applying the slice extraction for inputImage
    start[2] = maxSlice;

    MaskImageType::RegionType desiredRegion;
    desiredRegion.SetSize(size);
    desiredRegion.SetIndex(start);

    filter->SetExtractionRegion(desiredRegion);
    filter->SetInput(maskImage);
    filter->Update();
    mask2DImage = filter->GetOutput();


    // applying the slice extraction for inputImage
    {
        // InputImageType::RegionType inputRegion = inputImage->GetLargestPossibleRegion();
        // InputImageType::IndexType start = inputRegion.GetIndex();
        // InputImageType::SizeType size = inputRegion.GetSize();
        //
        // size[2] = 0;
        // start[2] = maxSlice;

        using FilterType = itk::ExtractImageFilter<InputImageType, InputImage2DType>;
        FilterType::Pointer filter = FilterType::New();
        // filter->InPlaceOn();
        filter->SetDirectionCollapseToSubmatrix();

        cout << start << size << endl;

        InputImageType::RegionType desiredRegion;
        desiredRegion.SetSize(size);
        desiredRegion.SetIndex(start);

        filter->SetExtractionRegion(desiredRegion);
        filter->SetInput(inputImage);
        filter->Update();
        input2DImage = filter->GetOutput();
    }

    return 0;
} // ExtractLargestAreaSlice

int
main(int argc, char * argv[])
{
    bool is2DShapeIntensity = false; // 2
    bool isIntensity        = false; // i
    bool isShapeIntensity   = false; // s
    bool isGeometry         = false; // g
    bool isGLCM = false;             // c
    bool isGLRM = false;             // r
    string featureSelect = "2isgcr";


    if (argc < 3) {
        cerr << "Missing Parameters " << endl;
        cerr << "Usage = " << argv[0];
        cerr << " InputImage LabelImage FeatureFile Label={1} FeatureSelect={2isgcr} NumberOfBins={64}" << endl;
        cerr << " 2 - 2DShapeIntensity, i - Intensity, s - Shape Intensity, g - Geometry, c - GLCM, r - GLRM" << endl;
        return EXIT_FAILURE;
    }

    LabelPixelType selectedLabelValue = maskValue;
    if (argc > 4) {
        selectedLabelValue = atoi(argv[4]);
    }

    if (argc > 5) {
        featureSelect = argv[5];
    }

    for (int i = 0; i < featureSelect.length(); i++) {
        switch (featureSelect[i]) {
            case '2':
                is2DShapeIntensity = true; // 2
                break;
            case 'i':
                isIntensity = true; // i
                break;
            case 's':
                isShapeIntensity = true; // s
                break;
            case 'g':
                isGeometry = true; // g
                break;
            case 'c':
                isIntensity = true; // c needs i
                isGLCM      = true; // c
                break;
            case 'r':
                isIntensity = true; // r needs i
                isGLRM      = true; // r
                break;
        }
    }

    string inputImageName  = argv[1];
    string labelImageName  = argv[2];
    string featureFileName = argv[3];

    if (argc > 6) {
        NumberOfBins = atoi(argv[6]);
    }


    cout << "Input Image Name = " << inputImageName << endl;
    cout << "Label Image Name = " << labelImageName << endl;
    cout << "Feature File Name = " << featureFileName << endl;


    /////////To Read CT images ////////////////////
    InputImageType::Pointer inputImage = ReadImageFile<InputImageType>(argv[1]);

    ///////// To Read the Mask images ////////////////////
    cout << "Get the mask from Input" << endl;
    LabelImageType::Pointer labelImage = ReadImageFile<LabelImageType>(argv[2]);

    // To check coordinate of the inputImage Image
    SpacingType inputImageSpacing     = inputImage->GetSpacing();
    OriginType inputImageOrigin       = inputImage->GetOrigin();
    RegionType inputImageRegion       = inputImage->GetLargestPossibleRegion();
    SizeType inputImageSize           = inputImageRegion.GetSize();
    DirectionType inputImageDirection = inputImage->GetDirection();

    // To check coordinate For mask Image
    SpacingType labelImageSpacing     = labelImage->GetSpacing();
    OriginType labelImageOrigin       = labelImage->GetOrigin();
    RegionType labelImageRegion       = labelImage->GetLargestPossibleRegion();
    SizeType labelImageRegionSize     = labelImageRegion.GetSize();
    DirectionType labelImageDirection = labelImage->GetDirection();

    /////////////////////////////

    cout << "Input Image Spacing = " << inputImageSpacing << endl;
    cout << "Input Image Origin = " << inputImageOrigin << endl;
    cout << "Input Image Size = " << inputImageSize << endl;
    cout << "Input Image Direction = " << inputImageDirection << endl << endl << endl;


    cout << "Label Image Spacing = " << labelImageSpacing << endl;
    cout << "Label Image Origin = " << labelImageOrigin << endl;
    cout << "Label Image RegionSize = " << labelImageRegionSize << endl;
    cout << "Label Image Direction = " << labelImageDirection << endl << endl << endl;

    // using TransformType = itk::IdentityTransform<double, Dimension>;

    /*
     * using ResampleImageFilterType = itk::ResampleImageFilter<InputImageType, InputImageType>;
     * ResampleImageFilterType::Pointer resampler = ResampleImageFilterType::New();
     * resampler->SetInput(inputImage);
     * resampler->SetSize(labelImageRegionSize);
     * resampler->SetOutputOrigin(labelImageOrigin);
     * resampler->SetOutputSpacing(labelImageSpacing);
     * resampler->SetOutputDirection(labelImageDirection);
     * //resampler->SetTransform(TransformType::New());
     * resampler->UpdateLargestPossibleRegion();
     * inputImage = resampler->GetOutput();
     *
     * inputImageSpacing = inputImage->GetSpacing();
     * inputImageOrigin  = inputImage->GetOrigin();
     * inputImageRegion  = inputImage->GetLargestPossibleRegion();
     * inputImageSize    = inputImageRegion.GetSize();
     *
     * cout << "Resampled Input Image Spacing = " << inputImageSpacing << endl;
     * cout << "Resampled Input Image Origin = " << inputImageOrigin << endl;
     * cout << "Resampled Input Image Size = " << inputImageSize << endl;
     * cout << "Resampled Input Image Direction = " << inputImageDirection << endl<< endl << endl;
     *
     * WriteImageFile<InputImageType>(inputImage, "test.nrrd");*/

    // Image ROI
    InputImageType::RegionType inputRegion;
    MaskImageType::RegionType maskRegion;

    MaskImageType::Pointer maskImage = MaskImageType::New();
    maskImage->SetRegions(labelImage->GetRequestedRegion());
    maskImage->CopyInformation(labelImage);
    maskImage->Allocate();

    maskImage = GetMaskImage(labelImage, selectedLabelValue);

    // For Island Removing
    MaskImageType::Pointer maskImageIslandRemoved;
    {
        // using LabelOpeningType = itk::LabelShapeKeepNObjectsImageFilter< MaskImageType >;
        using LabelOpeningType = itk::BinaryShapeKeepNObjectsImageFilter<MaskImageType>;

        LabelOpeningType::Pointer opening = LabelOpeningType::New();
        opening->SetInput(maskImage);
        opening->SetBackgroundValue(0);
        opening->SetFullyConnected(true);
        opening->SetForegroundValue(1);
        opening->SetNumberOfObjects(1);
        // opening->SetReverseOrdering( false );
        opening->SetAttribute(LabelOpeningType::LabelObjectType::NUMBER_OF_PIXELS); //  extract the largest object  NUMBER_OF_PIXELS
        opening->Update();

        maskImageIslandRemoved = opening->GetOutput();
    }
    maskImage = maskImageIslandRemoved;

    maskRegion = GetRoi(maskImage);
    // ExpandRoi(maskImage, maskRegion, 1);
    BoundingCheck(maskImage, inputImage, maskRegion, inputRegion);

    MaskImageType::IndexType roiStart;
    MaskImageType::IndexType roiEnd;

    RegionToIndex(maskRegion, roiStart, roiEnd);

    inputImage = ApplyRoi<InputImageType>(inputImage, inputRegion);
    maskImage  = ApplyRoi<MaskImageType>(maskImage, maskRegion);

    // To check coordinate For mask Image
    SpacingType maskImageSpacing     = maskImage->GetSpacing();
    OriginType maskImageOrigin       = maskImage->GetOrigin();
    RegionType maskImageRegion       = maskImage->GetLargestPossibleRegion();
    SizeType maskImageRegionSize     = maskRegion.GetSize();
    DirectionType maskImageDirection = maskImage->GetDirection();

    cout << "Mask Image Spacing = " << maskImageSpacing << endl;
    cout << "Mask Image Origin = " << maskImageOrigin << endl;
    cout << "Mask Image RegionSize = " << maskImageRegionSize << endl;
    cout << "Mask Image Direction = " << maskImageDirection << endl << endl << endl;

    inputImage->SetSpacing(maskImage->GetSpacing());
    inputImage->SetOrigin(maskImage->GetOrigin());
    inputImage->SetDirection(maskImage->GetDirection());

    InternalImageType::Pointer maskImageFloat;
    using CastToInputFilterType = itk::CastImageFilter<MaskImageType, InputImageType>;
    CastToInputFilterType::Pointer toInputFilter = CastToInputFilterType::New();
    toInputFilter->SetInput(maskImage);
    toInputFilter->Update();
    maskImageFloat = toInputFilter->GetOutput();


    // InternalImageType::Pointer inputImageFloat = InternalImageType::New();
    //    inputImageFloat->SetRegions( inputImage->GetRequestedRegion() );
    //    inputImageFloat->CopyInformation( inputImage );
    //    inputImageFloat->Allocate();

    //    {
    //        // Changing the data type of Input Image (unsigned int)  to float

    //        using inputImageIteratorFloatType = itk::ImageRegionIterator< InternalImageType>;
    //        inputImageIteratorFloatType inputImageFloatIt( inputImageFloat, inputImageFloat->GetRequestedRegion() );

    //        using inputImageIteratorType = itk::ImageRegionConstIterator< InputImageType>;
    //        inputImageIteratorType inputImageit( inputImage, inputImage->GetRequestedRegion() );

    //        inputImageit.GoToBegin();
    //        inputImageFloatIt.GoToBegin();

    //        while ( !inputImageit.IsAtEnd() )
    //        {
    //            inputImageFloatIt.Set(  inputImageit.Get() );
    //            ++inputImageit;
    //            ++inputImageFloatIt;
    //        }
    //    }


    // Output file open
    ofstream outfile((featureFileName).c_str(), ios::out);
    if (!outfile) {
        cerr << "Can't open output file " << featureFileName << endl;
        exit(1);
    }

    #if 0
    using SUVCastFilterType = itk::CastImageFilter<InputImageType, LabelImageType>;
    SUVCastFilterType::Pointer SUVCaster = SUVCastFilterType::New();

    SUVCaster->SetInput(inputImage);
    SUVCaster->Update();
    #endif


    #if 0
    {
        MaskWriterType::Pointer binaryWriter = MaskWriterType::New();

        binaryWriter->SetFileName("test0.dcm");
        binaryWriter->SetInput(maskImage);
        binaryWriter->Update();
    }
    #endif


    // maskImage = StentRemoval(inputImage, maskImage, 500);

    ///////////////////////////////////////////////////////////////////////////
    // Compute 2D shape and intensity feature uisng LabelMap   ////////////////////
    ////////////////////////////////////////////////////////////////////////////
    if (is2DShapeIntensity) {
        // extract the slice having the largest area
        MaskImage2DType::Pointer mask2DImage;
        InputImage2DType::Pointer input2DImage;
        ExtractLargestAreaSlice(inputImage, maskImage, input2DImage, mask2DImage);
        CalcShapeAndIntensityFeatures(input2DImage, mask2DImage, outfile);
    }


    /////////////////////////////////////////////////////////////////
    // Compute Intensity Features for three SUV images using itkLabelStatisticsImageFilter /////////
    ///////////////////////////////////////////////////////////////
    InputPixelType inputMin;
    InputPixelType inputMax;
    if (isIntensity)
        CalcIntensityFeatures(inputImage, maskImage, inputMin, inputMax, outfile);

    /////////////////////////////////////////////////////////////////////
    // Compute Geometry feature using LabelGeometryImageFilter  ////
    /////////////////////////////////////////////////////////////////////
    if (isGeometry)
        CalcGeometryFeatures(inputImage, maskImage, outfile);


    ///////////////////////////////////////////////////////////////////////////
    // Compute shape and intensity feature uisng LabelMap   ////////////////////
    ////////////////////////////////////////////////////////////////////////////
    if (isShapeIntensity)
        CalcShapeAndIntensityFeatures(inputImage, maskImage, outfile);


    ////////////////////////////////////////////////////////////////////
    ///// Compute Feature for the Input images using ScalarImageToTextureFeaturesFilter /////
    ////////////////////////////////////////////////////////////////////
    if (isIntensity && isGLCM)
        CalcGlcmFeatures(inputImage, maskImageFloat, inputMin, inputMax, outfile);


    ////////////////////////////////////////////////////////////////////
    ///// Compute Feature for the Input images using ScalarImageToRunLengthFeaturesFilter /////
    ////////////////////////////////////////////////////////////////////
    if (isIntensity && isGLRM) {
        MaskImageType::PointType pointMin;
        MaskImageType::PointType pointMax;

        maskImage->TransformIndexToPhysicalPoint(roiStart, pointMin);
        maskImage->TransformIndexToPhysicalPoint(roiEnd, pointMax);

        double distMax = pointMin.EuclideanDistanceTo(pointMax);

        CalcGlrmFeatures(inputImage, maskImageFloat, inputMin, inputMax, distMax, outfile);
    }

    #if 0
    binaryWriter->SetFileName("test1.dcm");
    binaryWriter->SetInput(maskImage);
    binaryWriter->Update();
    #endif

    ////////////////////////////////////////////////////////////////////
    ///// Compute Glucosity Feature /////
    ////////////////////////////////////////////////////////////////////
    #if 0
    {
        outfile << "Glucosity = " << (labelObjectS->GetPhysicalSize()) * (labelObjectS->GetMean()) << endl;
    }
    #endif

    outfile.close();

    return EXIT_SUCCESS;
} // main


<|start_filename|>Tools/ITKUtils/itkMaskedSpatialObjectToImageFilter.hxx<|end_filename|>
/*=========================================================================

Program:   Insight Segmentation & Registration Toolkit
Module:    $RCSfile: itkMaskedSpatialObjectToImageFilter.cxx,v $
Language:  C++
Date:      $Date: 2006/07/19 20:06:17 $
Version:   $Revision: 1.00 $

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

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.

=========================================================================*/
#ifndef __itkMaskedSpatialObjectToImageFilter_cxx
#define __itkMaskedSpatialObjectToImageFilter_cxx

#include "itkMaskedSpatialObjectToImageFilter.h"
#include "itkSpatialObjectToImageFilter.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkNumericTraits.h"
#include "itkImageRegionConstIterator.h"
#include "itkProgressReporter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkBinaryDilateImageFilter.h"
#include "itkBinaryBallStructuringElement.h"
#include "itkResampleImageFilter.h"
#include "itkCenteredAffineTransform.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkBinaryThresholdImageFilter.h"

namespace itk
{

/** Constructor */
template <class TInputSpatialObject, class TOutputImage>
MaskedSpatialObjectToImageFilter<TInputSpatialObject, TOutputImage>::MaskedSpatialObjectToImageFilter()
{
  this->SetNumberOfRequiredInputs(1);
  this->m_ChildrenDepth = 1;
  this->m_Size.Fill(0);

  for (unsigned int i = 0; i < itkGetStaticConstMacro(OutputImageDimension); i++)
  {
    this->m_Spacing[i] = 1.0;
    this->m_Origin[i] = 0.0;
  }

  this->m_InsideValue = 0;
  this->m_OutsideValue = 0;
  this->m_UseObjectValue = false;

  this->m_UserSuppliedMask = false;
  this->m_MaskResampleFactor = 4.0;
  this->m_MaskDilationSize = 2;
}

/** Destructor */
template <class TInputSpatialObject, class TOutputImage>
MaskedSpatialObjectToImageFilter<TInputSpatialObject, TOutputImage>::~MaskedSpatialObjectToImageFilter()
{
}

/** Set the Input SpatialObject */
template <class TInputSpatialObject, class TOutputImage>
void MaskedSpatialObjectToImageFilter<TInputSpatialObject, TOutputImage>::SetInput(const InputSpatialObjectType *input)
{
  // Process object is not const-correct so the const_cast is required here
  this->ProcessObject::SetNthInput(0,
                                   const_cast<InputSpatialObjectType *>(input));
}

/** Set the Input Mask image */
template <class TInputSpatialObject, class TOutputImage>
void MaskedSpatialObjectToImageFilter<TInputSpatialObject, TOutputImage>::SetMask(const MaskImageType *mask)
{
  if (mask)
  {
    this->m_UserSuppliedMask = true;

    // Const-correct the mask image
    this->ProcessObject::SetNthInput(1,
                                     const_cast<MaskImageType *>(mask));
  }
}

/** Get the input Spatial Object */
template <class TInputSpatialObject, class TOutputImage>
const typename MaskedSpatialObjectToImageFilter<TInputSpatialObject, TOutputImage>::InputSpatialObjectType *
MaskedSpatialObjectToImageFilter<TInputSpatialObject, TOutputImage>::GetInput(void)
{
  if (this->GetNumberOfInputs() < 1)
  {
    return 0;
  }

  return static_cast<const TInputSpatialObject *>(this->ProcessObject::GetInput(0));
}

/** Get the input Mask image */
template <class TInputSpatialObject, class TOutputImage>
const typename MaskedSpatialObjectToImageFilter<TInputSpatialObject, TOutputImage>::MaskImageType *
MaskedSpatialObjectToImageFilter<TInputSpatialObject, TOutputImage>::GetMask(void)
{
  if (this->GetNumberOfInputs() < 2)
  {
    return 0;
  }

  return static_cast<const MaskImageType *>(this->ProcessObject::GetInput(1));
}

/** GenerateData */
template <class TInputSpatialObject, class TOutputImage>
void MaskedSpatialObjectToImageFilter<TInputSpatialObject, TOutputImage>::GenerateData(void)
{
  this->SetupOutputImage();

  if (!this->m_UserSuppliedMask)
    //We must generate the mask ourselves
    this->GenerateMask();

  this->GenerateDataUsingMask();
}

/** SetupOutputImage */
template <class TInputSpatialObject, class TOutputImage>
void MaskedSpatialObjectToImageFilter<TInputSpatialObject, TOutputImage>::SetupOutputImage(void)
{
  unsigned int i;

  // Get the input, and output pointers
  const InputSpatialObjectType *InputObject = this->GetInput();
  OutputImagePointer OutputImage = this->GetOutput();

  // Generate the image
  SizeType size;

#if ITK_VERSION_MAJOR == 5
  InputObject->ComputeFamilyBoundingBox(this->m_ChildrenDepth);
  for (i = 0; i < this->ObjectDimension; i++)
  {
    size[i] = static_cast<SizeValueType>(
        InputObject->GetFamilyBoundingBoxInWorldSpace()->GetMaximum()[i] - InputObject->GetFamilyBoundingBoxInWorldSpace()->GetMinimum()[i]);
  }
#else
  InputObject->ComputeBoundingBox();
  for (unsigned int i = 0; i < itkGetStaticConstMacro(ObjectDimension); i++)
  {
    size[i] = (long unsigned int)(InputObject->GetBoundingBox()->GetMaximum()[i] - InputObject->GetBoundingBox()->GetMinimum()[i]);
  }
#endif

  typename OutputImageType::IndexType index;
  index.Fill(0);
  typename OutputImageType::RegionType region;

  // If the size of the output has been explicitly specified, the filter
  // will set the output size to the explicit size, otherwise the size from the
  // spatial
  // object's bounding box will be used as default.

  bool specified = false;
  for (i = 0; i < this->OutputImageDimension; i++)
  {
    if (this->m_Size[i] != 0)
    {
      specified = true;
      break;
    }
  }

  if (specified)
  {
    region.SetSize(this->m_Size);
  }
  else
  {
    region.SetSize(size);
  }
  region.SetIndex(index);

  OutputImage->SetLargestPossibleRegion(region); //
  OutputImage->SetBufferedRegion(region);        // set the region
  OutputImage->SetRequestedRegion(region);       //
  OutputImage->SetSpacing(this->m_Spacing);      // set spacing
  OutputImage->SetOrigin(this->m_Origin);        //   and origin
  OutputImage->SetDirection(this->m_Direction);
  OutputImage->Allocate(); // allocate the image
}

/** Generate Mask */
template <class TInputSpatialObject, class TOutputImage>
void MaskedSpatialObjectToImageFilter<TInputSpatialObject, TOutputImage>::GenerateMask(void)
{
  // NOTE: The mask is store in Inputs(1)
  MaskImagePointer maskImage = nullptr;

  //Convert SpatialObject to low-resolution mask
  //NOTE: We set the InsideValue and Outside values to max() and zero, and
  //      tell the filter to not use the object value to ensure that
  //      grey-scale SpatialObjects are converted to a binary mask.
  using SpatialObjectToImageFilterType = itk::SpatialObjectToImageFilter<TInputSpatialObject, MaskImageType>;
  typename SpatialObjectToImageFilterType::Pointer filterConvertMask =
      SpatialObjectToImageFilterType::New();
  filterConvertMask->SetInput(this->GetInput());
  typename SpatialObjectToImageFilterType::SizeType maskSize;       //Set Size
  typename SpatialObjectToImageFilterType::SpacingType maskSpacing; //Set Spacing
  for (unsigned int dim = 0; dim < itkGetStaticConstMacro(OutputImageDimension); dim++)
  {
    maskSize[dim] = (unsigned int)(this->m_Size[dim] / this->m_MaskResampleFactor);
    maskSpacing[dim] = this->m_Spacing[dim] * ((double)this->m_Size[dim] / (double)maskSize[dim]);
  }
  //std::cout << this->m_Size << " " << maskSize << std::endl;
  //std::cout << (SpacingType)this->m_Spacing << " " << maskSpacing << std::endl;
  filterConvertMask->SetOrigin(this->m_Origin);
  filterConvertMask->SetSize(maskSize);
  filterConvertMask->SetSpacing(maskSpacing);
  filterConvertMask->SetInsideValue(itk::NumericTraits<MaskPixelType>::max());
  filterConvertMask->SetOutsideValue(itk::NumericTraits<MaskPixelType>::Zero);
  filterConvertMask->SetUseObjectValue(false);
  filterConvertMask->Update();

  //Create structuring element
  using BinaryBallElementType = itk::BinaryBallStructuringElement<MaskPixelType, itkGetStaticConstMacro(OutputImageDimension)>;
  BinaryBallElementType kernel;
  kernel.SetRadius(this->m_MaskDilationSize);
  kernel.CreateStructuringElement();

  //Binary dilate the resampled Mask image
  using BinaryDilateImageFilterType = itk::BinaryDilateImageFilter<MaskImageType, MaskImageType, BinaryBallElementType>;
  typename BinaryDilateImageFilterType::Pointer filterDilateMask = BinaryDilateImageFilterType::New();
  filterDilateMask->SetInput(filterConvertMask->GetOutput());
  filterDilateMask->SetKernel(kernel);
  filterDilateMask->SetDilateValue(itk::NumericTraits<MaskPixelType>::max());
  filterDilateMask->Update();

  //Resample mask to full size
  using ResampleFilterType = itk::ResampleImageFilter<MaskImageType, MaskImageType, MaskCoordRep>;
  using TransformType = itk::CenteredAffineTransform<MaskCoordRep, itkGetStaticConstMacro(OutputImageDimension)>;
  using InterpolatorType = itk::NearestNeighborInterpolateImageFunction<MaskImageType, MaskCoordRep>;

  typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
  typename TransformType::Pointer transform = TransformType::New();
  transform->SetIdentity();
  typename ResampleFilterType::Pointer filterMaskResample = ResampleFilterType::New();
  filterMaskResample->SetTransform(transform);
  filterMaskResample->SetInterpolator(interpolator);
  filterMaskResample->SetInput(filterDilateMask->GetOutput());
  double finalMaskOrigin[itkGetStaticConstMacro(OutputImageDimension)]; //Set Origin
  typename ResampleFilterType::SizeType finalMaskSize;                  //Set Size
  typename ResampleFilterType::SpacingType finalMaskSpacing;            //Set Spacing
  for (unsigned int dim = 0; dim < itkGetStaticConstMacro(OutputImageDimension); dim++)
  {
    finalMaskOrigin[dim] = this->m_Origin[dim];
    finalMaskSize[dim] = this->m_Size[dim];
    finalMaskSpacing[dim] = this->m_Spacing[dim];
  }
  //std::cout << filterDilateMask->GetOutput()->GetLargestPossibleRegion().GetSize() << " " << finalMaskSize << std::endl;
  //std::cout << filterDilateMask->GetOutput()->GetSpacing() << " " << finalMaskSpacing << std::endl;

  filterMaskResample->SetOutputOrigin(finalMaskOrigin);
  filterMaskResample->SetSize(finalMaskSize);
  filterMaskResample->SetOutputSpacing(finalMaskSpacing);
  filterMaskResample->Update();

  //Set mask image
  maskImage = filterMaskResample->GetOutput();
  this->SetMask(maskImage);
}

/** GenerateData (using Mask) */
template <class TInputSpatialObject, class TOutputImage>
void MaskedSpatialObjectToImageFilter<TInputSpatialObject, TOutputImage>::GenerateDataUsingMask(void)
{
  itkDebugMacro(<< "MaskedSpatialObjectToImageFilter::Update() called");

  //Get the input, mask and output pointers
  const InputSpatialObjectType *InputObject = this->GetInput();
  OutputImagePointer OutputImage = this->GetOutput();
  const MaskImageType *MaskImage = this->GetMask();

  //Get output region
  typename OutputImageType::RegionType region =
      OutputImage->GetLargestPossibleRegion();

  //Setup Progress reporter
  ProgressReporter progress(this, 0, region.GetNumberOfPixels());

  using OutputIteratorType = itk::ImageRegionIteratorWithIndex<OutputImageType>;
  using MaskIteratorType = itk::ImageRegionConstIterator<MaskImageType>;

  OutputIteratorType itOutput(OutputImage, region);
  MaskIteratorType itMask(MaskImage, region);

  itk::Point<double, itkGetStaticConstMacro(ObjectDimension)> point;

  while (!itOutput.IsAtEnd())
  {
    //Check if the mask is "ON" for this point
    if (itMask.Get() > NumericTraits<MaskPixelType>::Zero)
    {
      // ValueAt requires the point to be in physical coordinate i.e
      for (unsigned int i = 0; i < itkGetStaticConstMacro(ObjectDimension); i++)
      {
        point[i] = (itOutput.GetIndex()[i] * this->m_Spacing[i]) + this->m_Origin[i];
      }
      double val = 0.0;
#if ITK_VERSION_MAJOR == 5
      InputObject->ValueAtInObjectSpace(point, val, TInputSpatialObject::MaximumDepth);
#else
      InputObject->ValueAt(point, val, TInputSpatialObject::MaximumDepth);
#endif

      if (this->m_InsideValue != 0 || this->m_OutsideValue != 0)
      {
        if (val)
        {
          if (this->m_UseObjectValue)
          {
            itOutput.Set(static_cast<ValueType>(val));
          }
          else
          {
            itOutput.Set(this->m_InsideValue);
          }
        }
        else
        {
          itOutput.Set(this->m_OutsideValue);
        }
      }
      else
      {
        itOutput.Set(static_cast<ValueType>(val));
      }
    }
    else
    {
      //The mask is "OFF" for this region - set as outside
      itOutput.Set(this->m_OutsideValue);
    }
    ++itOutput;
    ++itMask;
    progress.CompletedPixel();
  }

  itkDebugMacro(<< "MaskedSpatialObjectToImageFilter::Update() finished");

} // end update function

template <class TInputSpatialObject, class TOutputImage>
void MaskedSpatialObjectToImageFilter<TInputSpatialObject, TOutputImage>::PrintSelf(std::ostream &os, Indent indent) const
{
  Superclass::PrintSelf(os, indent);
  os << indent << "Size : " << this->m_Size << std::endl;
  os << indent << "Spacing : " << (SpacingType)this->m_Spacing << std::endl;
  os << indent << "Children depth : " << this->m_ChildrenDepth << std::endl;
  os << indent << "Inside Value : " << this->m_InsideValue << std::endl;
  os << indent << "Outside Value : " << this->m_OutsideValue << std::endl;
  if (this->m_UseObjectValue)
  {
    os << indent << "Using Object Value : ON" << std::endl;
  }
  else
  {
    os << indent << "Using Object Value : OFF" << std::endl;
  }
  if (this->m_UserSuppliedMask)
  {
    os << indent << "User Supplied Mask : YES" << std::endl;
  }
  else
  {
    os << indent << "User Supplied Mask : NO" << std::endl;
  }
}

} // end namespace itk

#endif


<|start_filename|>Tools/ContourTools/HoleGenerator.cxx<|end_filename|>
#include "ITKUtils.h"

#include "itkSpatialObjectToImageFilter.h"
#include "itkEllipseSpatialObject.h"

#include "itkImageDuplicator.h"

using namespace std;

int main( int argc, char *argv[] )
{
  if ( argc < 7 )
  {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " InputLabelImageFile SeedPoint[x y z] HoleSize OutputPrefix" << std::endl;
    return EXIT_FAILURE;
  }

  const unsigned int x = atoi( argv[2] );
  const unsigned int y = atoi( argv[3] );
  const int sliceNumber = atoi( argv[4] );
  const double holeD = atof( argv[5] );
  const double holeR = holeD / 2;
  const string outputPrefix = argv[6];


  /////////To Read Label images ////////////////////
  LabelImageType::Pointer labelImage = ReadImageFile<MaskImageType>(argv[1]);
  LabelImageType::Pointer sphereImage = LabelImageType::New();
  LabelImageType::Pointer holedImage;

  sphereImage->CopyInformation(labelImage);
  sphereImage->SetBufferedRegion(labelImage->GetBufferedRegion());
  sphereImage->Allocate();
  sphereImage->FillBuffer( 0 );

  using DuplicatorType = itk::ImageDuplicator< LabelImageType >;
  DuplicatorType::Pointer duplicator = DuplicatorType::New();
  duplicator->SetInputImage(labelImage);
  duplicator->Update();
  holedImage = duplicator->GetOutput();;

  // set ROI
  LabelImageType::PointType origin = labelImage->GetOrigin();
  LabelImageType::IndexType index = labelImage->GetLargestPossibleRegion().GetIndex();
  LabelImageType::SizeType size = labelImage->GetLargestPossibleRegion().GetSize();
  LabelImageType::SpacingType spacing = labelImage->GetSpacing();
  LabelImageType::SizeType roiSize;
  LabelImageType::IndexType roiStart;
  LabelImageType::IndexType roiEnd;
  LabelImageType::IndexType center;

  roiSize[0] = static_cast<unsigned int>((holeD + spacing[2]*2) / spacing[0]);
  roiSize[1] = static_cast<unsigned int>((holeD + spacing[2]*2) / spacing[1]);
  roiSize[2] = static_cast<unsigned int>((holeD + spacing[2]*2) / spacing[2]);

  center[0] = x;
  center[1] = y;

  if (sliceNumber > 0)
    center[2] = size[2] - sliceNumber - 1;
  else
    center[2] = -sliceNumber;

  std::cout << center << std::endl;

  roiStart[0] = center[0] - roiSize[0] / 2;
  roiStart[1] = center[1] - roiSize[1] / 2;
  roiStart[2] = center[2] - roiSize[2] / 2;
  roiEnd[0] = roiStart[0] + roiSize[0];
  roiEnd[1] = roiStart[1] + roiSize[1];
  roiEnd[2] = roiStart[2] + roiSize[2];

  std::cout << roiSize << std::endl;
  std::cout << roiStart << std::endl;
  std::cout << roiEnd << std::endl;

  // make a sphere
  {
    using EllipseType = itk::EllipseSpatialObject<3>;
    using SpatialObjectToImageFilterType = itk::SpatialObjectToImageFilter<EllipseType, LabelImageType>;

    using TransformType = EllipseType::TransformType;

    SpatialObjectToImageFilterType::Pointer sphereFilter = SpatialObjectToImageFilterType::New();

    TransformType::Pointer transform = TransformType::New();

    transform->SetIdentity();

    TransformType::OutputVectorType  translation;
    translation[0] = roiSize[0] * spacing[0] / 2.0;
    translation[1] = roiSize[1] * spacing[1] / 2.0;
    translation[2] = roiSize[2] * spacing[2] / 2.0;
    transform->Translate( translation, false );

    std::cout << translation << std::endl;

    EllipseType::ArrayType radiusArray;
    EllipseType::Pointer sphere = EllipseType::New();
    radiusArray[0] = holeR;
    radiusArray[1] = holeR;
    radiusArray[2] = holeR;

    sphere->SetRadius(radiusArray);
    sphere->SetDefaultInsideValue(maskValue);
    sphere->SetDefaultOutsideValue(0);
    sphere->SetObjectToParentTransform(transform);

    sphereFilter->SetSize(roiSize);
    sphereFilter->SetSpacing(spacing);
    sphereFilter->SetInput(sphere);
    sphereFilter->SetUseObjectValue(true);
    sphereFilter->SetOutsideValue(0);

    try
    {
      sphereFilter->Update();

      itk::ImageRegionIteratorWithIndex<LabelImageType> sphere(sphereImage, sphereImage->GetBufferedRegion() );
      itk::ImageRegionIterator<LabelImageType> holed(holedImage, holedImage->GetBufferedRegion() );
      itk::ImageRegionConstIterator<LabelImageType> sphereROI(sphereFilter->GetOutput(), sphereFilter->GetOutput()->GetBufferedRegion());

      for (sphere.GoToBegin(), holed.GoToBegin(), sphereROI.GoToBegin(); !sphere.IsAtEnd();
           ++sphere,++holed)
      {
        LabelImageType::IndexType idx = sphere.GetIndex();

        if (idx[0] >= roiStart[0] && idx[1] >= roiStart[1] && idx[2] >= roiStart[2] &&
            idx[0] < roiEnd[0] && idx[1] < roiEnd[1] && idx[2] < roiEnd[2])
        {
          LabelImageType::PixelType color = sphereROI.Get();

          ++sphereROI;

          if (color > 0)
          {
            sphere.Set(maskValue); // sphere
            holed.Set(0); // hole
          }
        }
      }
    }
    catch ( itk::ExceptionObject &excp )
    {
      std::cerr << excp << std::endl;
      return EXIT_FAILURE;
    }
  }


  try
  {
    WriteImageFile<LabelImageType>(sphereImage, outputPrefix + "_sphere_" + argv[5] + "-label.nrrd");
    WriteImageFile<LabelImageType>(holedImage, outputPrefix + "_holed_" + argv[5] + "-label.nrrd");
  }
  catch (itk::ExceptionObject &excp)
  {
    std::cerr << "Exception thrown while writing the series " << std::endl;
    std::cerr << excp << std::endl;
    return EXIT_FAILURE;
  }


  std::cout << "Saved the label image" << std::endl;


  return EXIT_SUCCESS;
}


<|start_filename|>install_modules.cmd<|end_filename|>
conda config --set always_yes yes --set changeps1 no
conda update -q conda
conda install --yes python=$TRAVIS_PYTHON_VERSION pandas scipy numpy ipython matplotlib
conda install --yes -c bioconda ruffus
conda install --yes -c SimpleITK SimpleITK 
REM conda install pyqt pyside


<|start_filename|>Tools/ITKUtils/ITKUtils.cxx<|end_filename|>
#if defined(_MSC_VER)
# pragma warning ( disable : 4786 )
#endif

#include "ITKUtils.h"

MaskImageType::Pointer
GetMaskImage(LabelImageType::Pointer labelImage)
{
    MaskImageType::Pointer maskImage = MaskImageType::New();
    maskImage->SetRegions(labelImage->GetRequestedRegion());
    maskImage->CopyInformation(labelImage);
    maskImage->Allocate();

    itk::ImageRegionIterator<LabelImageType> label(labelImage, labelImage->GetBufferedRegion());
    itk::ImageRegionIterator<MaskImageType> mask(maskImage, maskImage->GetBufferedRegion());

    for (label.GoToBegin(), mask.GoToBegin(); !label.IsAtEnd(); ++label, ++mask) {
        LabelImageType::PixelType color = label.Get();
        if (color > 0) {
            mask.Set(maskValue);
        } else {
            mask.Set(0);
        }
    }

    return maskImage;
}

MaskImageType::Pointer
GetMaskImage(LabelImageType::Pointer labelImage, LabelPixelType selectedLabelValue)
{
    MaskImageType::Pointer maskImage = MaskImageType::New();
    maskImage->SetRegions(labelImage->GetRequestedRegion());
    maskImage->CopyInformation(labelImage);
    maskImage->Allocate();

    itk::ImageRegionIterator<LabelImageType> label(labelImage, labelImage->GetBufferedRegion());
    itk::ImageRegionIterator<MaskImageType> mask(maskImage, maskImage->GetBufferedRegion());

    for (label.GoToBegin(), mask.GoToBegin(); !label.IsAtEnd(); ++label, ++mask) {
        LabelImageType::PixelType color = label.Get();
        if (color == selectedLabelValue) {
            mask.Set(maskValue);
        } else {
            mask.Set(0);
        }
    }

    return maskImage;
}

MaskImageType::RegionType
RoiIndexToRegion(MaskImageType::IndexType roiStart, MaskImageType::IndexType roiEnd)
{
    MaskImageType::RegionType outputReigion;
    MaskImageType::IndexType outputStart;
    MaskImageType::SizeType outputSize;

    for (unsigned n = 0; n < Dimension; n++) {
        outputStart[n] = roiStart[n];
        outputSize[n]  = roiEnd[n] - roiStart[n] + 1;
    }

    outputReigion.SetIndex(outputStart);
    outputReigion.SetSize(outputSize);

    return outputReigion;
}

void
RegionToIndex(MaskImageType::RegionType outputReigion, MaskImageType::IndexType &roiStart,
  MaskImageType::IndexType &roiEnd)
{
    MaskImageType::IndexType outputStart = outputReigion.GetIndex();
    MaskImageType::SizeType outputSize   = outputReigion.GetSize();

    for (unsigned n = 0; n < Dimension; n++) {
        roiStart[n] = outputStart[n];
        roiEnd[n]   = outputStart[n] + outputSize[n] + 1;
    }
}

MaskImageType::RegionType
ReadROI(string roiName)
{
    MaskImageType::RegionType roiRegion;
    MaskImageType::IndexType roiStart;
    MaskImageType::IndexType roiEnd;
    MaskImageType::SizeType roiSize;

    ifstream roiFile;
    char buf[256];
    roiFile.open(roiName.c_str(), ios::in);

    roiFile >> roiStart[0];
    roiFile >> roiStart[1];
    roiFile >> roiStart[2];
    roiFile.getline(buf, 256);
    cout << roiStart << endl;

    roiFile >> roiEnd[0];
    roiFile >> roiEnd[1];
    roiFile >> roiEnd[2];
    roiFile.getline(buf, 256);
    cout << roiEnd << endl;

    roiFile >> roiSize[0];
    roiFile >> roiSize[1];
    roiFile >> roiSize[2];
    roiFile.getline(buf, 256);
    cout << roiSize << endl;

    roiRegion.SetIndex(roiStart);
    roiRegion.SetSize(roiSize);

    roiFile.close();

    return roiRegion;
} // ReadROI

void
WriteROI(MaskImageType::RegionType inputRegion, string roiName)
{
    MaskImageType::IndexType roiStart;
    MaskImageType::IndexType roiEnd;
    RegionToIndex(inputRegion, roiStart, roiEnd);
    ofstream roiFile;

    roiFile.open(roiName.c_str(), ios::out);
    for (unsigned i = 0; i < Dimension; i++)
        roiFile << roiStart[i] << " ";
    roiFile << endl;
    for (unsigned i = 0; i < Dimension; i++)
        roiFile << roiEnd[i] << " ";
    roiFile << endl;
    for (unsigned i = 0; i < Dimension; i++)
        roiFile << inputRegion.GetSize()[i] << " ";
    roiFile << endl;
    roiFile.close();
}

MaskImageType::RegionType
GetRoi(MaskImageType::Pointer maskImage)
{
    MaskImageType::RegionType roiRegion;
    MaskImageType::IndexType roiStart;
    MaskImageType::IndexType roiEnd;

    roiStart[0] = 10000;
    roiStart[1] = 10000;
    roiStart[2] = 10000;
    roiEnd[0]   = 0;
    roiEnd[1]   = 0;
    roiEnd[2]   = 0;

    itk::ImageRegionIteratorWithIndex<MaskImageType> mask(maskImage, maskImage->GetBufferedRegion());
    for (mask.GoToBegin(); !mask.IsAtEnd(); ++mask) {
        MaskImageType::PixelType color = mask.Get();
        if (color > 0) {
            MaskImageType::IndexType idx = mask.GetIndex();
            for (unsigned i = 0; i < Dimension; i++) {
                if (idx[i] < roiStart[i]) {
                    roiStart[i] = idx[i];
                }
                if (idx[i] > roiEnd[i]) {
                    roiEnd[i] = idx[i];
                }
            }
        }
    }

    cout << "roiStart " << roiStart << ", roiEnd " << roiEnd << "" << endl;

    roiRegion = RoiIndexToRegion(roiStart, roiEnd);

    return roiRegion;
} // GetRoi

void
ExpandRoi(MaskImageType::Pointer maskImage, MaskImageType::RegionType &maskRegion, const double objectSize)
{
    RegionType maskImageRegion = maskImage->GetLargestPossibleRegion();
    SizeType maskImageSize     = maskImageRegion.GetSize();

    MaskImageType::IndexType maskStart = maskRegion.GetIndex();
    MaskImageType::SizeType maskSize   = maskRegion.GetSize();

    MaskImageType::IndexType index     = maskImage->GetBufferedRegion().GetIndex();
    MaskImageType::SizeType size       = maskImage->GetBufferedRegion().GetSize();
    MaskImageType::SpacingType spacing = maskImage->GetSpacing();

    cout << "maskStart " << maskStart << ", maskEnd " << maskStart + maskSize << "" << endl;

    for (unsigned i = 0; i < Dimension; i++) {
        MaskImageType::PixelType radius = static_cast<MaskImageType::PixelType>(objectSize / spacing[i]);

        maskStart[i] -= radius;
        maskSize[i]  += 2 * radius;

        maskStart[i] = (maskStart[i] < 0) ? 0 : maskStart[i];
        maskSize[i]  = (static_cast<unsigned int>(maskStart[i] + maskSize[i] - 1) < maskImageSize[i]) ?
          maskSize[i] : (maskImageSize[i] - maskStart[i]);
    }

    cout << "maskStart " << maskStart << ", maskEnd " << maskStart + maskSize << "" << endl;

    maskRegion.SetIndex(maskStart);
    maskRegion.SetSize(maskSize);
}

void
BoundingCheck(MaskImageType::Pointer maskImage, InputImageType::Pointer inputImage,
  MaskImageType::RegionType &maskRegion, InputImageType::RegionType &inputRegion)
{
    RegionType inputImageRegion = inputImage->GetLargestPossibleRegion();
    SizeType inputImageSize     = inputImageRegion.GetSize();

    MaskImageType::IndexType maskStart = maskRegion.GetIndex();
    MaskImageType::SizeType maskSize   = maskRegion.GetSize();
    MaskImageType::IndexType maskCenter;

    InputImageType::IndexType istart;
    InputImageType::SizeType isize;

    MaskImageType::IndexType ostart;
    MaskImageType::SizeType osize;

    MaskImageType::PointType ptRoiStart;
    maskImage->TransformIndexToPhysicalPoint(maskStart, ptRoiStart);
    inputImage->TransformPhysicalPointToIndex(ptRoiStart, istart);

    for (unsigned n = 0; n < Dimension; n++) {
        isize[n]  = maskSize[n];
        ostart[n] = maskStart[n];

        // Input image
        if (istart[n] < 0) {
            ostart[n] = ostart[n] - istart[n];
            isize[n]  = isize[n] + istart[n];
            istart[n] = 0;
        }
        isize[n] = (static_cast<unsigned int>(isize[n] + istart[n]) <= inputImageSize[n]) ?
          isize[n] : (inputImageSize[n] - istart[n] - 1);
        // Mask
        osize[n] = isize[n];

        maskCenter[n] = static_cast<MaskImageType::IndexValueType>(round(isize[n] / 2));
    }

    cout << "maskStart " << maskStart << ", maskEnd " << maskStart + maskSize << "" << endl;
    cout << "istart " << istart << ", isize " << isize << "" << endl;
    cout << "ostart " << ostart << ", osize " << osize << "" << endl;
    cout << "maskCenter " << maskCenter << endl;

    inputRegion.SetIndex(istart);
    inputRegion.SetSize(isize);

    maskRegion.SetIndex(ostart);
    maskRegion.SetSize(osize);
} // BoundingCheck

void
BoundingCheck(MaskImageType::Pointer maskImage, LabelImageType::Pointer inputImage,
  MaskImageType::RegionType &maskRegion, InputImageType::RegionType &inputRegion)
{
    RegionType inputImageRegion = inputImage->GetLargestPossibleRegion();
    SizeType inputImageSize     = inputImageRegion.GetSize();

    MaskImageType::IndexType maskStart = maskRegion.GetIndex();
    MaskImageType::SizeType maskSize   = maskRegion.GetSize();
    MaskImageType::IndexType maskCenter;

    LabelImageType::IndexType istart;
    LabelImageType::SizeType isize;

    MaskImageType::IndexType ostart;
    MaskImageType::SizeType osize;

    MaskImageType::PointType ptRoiStart;
    maskImage->TransformIndexToPhysicalPoint(maskStart, ptRoiStart);
    inputImage->TransformPhysicalPointToIndex(ptRoiStart, istart);

    for (unsigned n = 0; n < Dimension; n++) {
        isize[n]  = maskSize[n];
        ostart[n] = maskStart[n];

        // Input image
        if (istart[n] < 0) {
            ostart[n] = ostart[n] - istart[n];
            isize[n]  = isize[n] + istart[n];
            istart[n] = 0;
        }
        isize[n] = (static_cast<unsigned int>(isize[n] + istart[n]) <= inputImageSize[n]) ?
          isize[n] : (inputImageSize[n] - istart[n] - 1);
        // Mask
        osize[n] = isize[n];

        maskCenter[n] = static_cast<MaskImageType::IndexValueType>(round(isize[n] / 2));
    }

    cout << "maskStart " << maskStart << ", maskEnd " << maskStart + maskSize << "" << endl;
    cout << "istart " << istart << ", isize " << isize << "" << endl;
    cout << "ostart " << ostart << ", osize " << osize << "" << endl;
    cout << "maskCenter " << maskCenter << endl;

    inputRegion.SetIndex(istart);
    inputRegion.SetSize(isize);

    maskRegion.SetIndex(ostart);
    maskRegion.SetSize(osize);
} // BoundingCheck

MaskImageType::Pointer
unionImages(vector<MaskImageType::Pointer> inputMaskImages)
{
    MaskImageType::Pointer outputMaskImage = MaskImageType::New();
    outputMaskImage->SetRegions(inputMaskImages[0]->GetLargestPossibleRegion());
    outputMaskImage->CopyInformation(inputMaskImages[0]);
    outputMaskImage->Allocate();
    outputMaskImage->FillBuffer(0);

    int volume = 0;
    for (int i = 0; i < inputMaskImages.size(); i++) {
        itk::ImageRegionIterator<MaskImageType> in(inputMaskImages[i], inputMaskImages[i]->GetBufferedRegion());
        itk::ImageRegionIterator<MaskImageType> out(outputMaskImage, outputMaskImage->GetBufferedRegion());

        for (in.GoToBegin(), out.GoToBegin(); !in.IsAtEnd(); ++in, ++out) {
            MaskImageType::PixelType mask = in.Get();
            if (mask > 0) {
                out.Set(maskValue);
            }
        }
    }

    return outputMaskImage;
}


<|start_filename|>Tools/ContourTools/ROIGenerator.cxx<|end_filename|>
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkResampleImageFilter.h>
#include <itkCastImageFilter.h>
#include <itkCommand.h>

#include <itkRegionOfInterestImageFilter.h>
#include <itkExtractImageFilter.h>

// Operations on Images

#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIterator.h>
#include <vnl/vnl_math.h>

// For Feature
#include <itkRescaleIntensityImageFilter.h>


#include <itkMinimumMaximumImageCalculator.h>
#include <itkStatisticsImageFilter.h>

#include <itkScalarImageToTextureFeaturesFilter.h>
#include <itkScalarImageToRunLengthFeaturesFilter.h>

#include <itkLabelGeometryImageFilter.h>
#include <itkConnectedComponentImageFilter.h>
#include <itkRelabelComponentImageFilter.h>
#include <itkLabelStatisticsImageFilter.h>

// For Label Object Representation
#include <itkLabelObject.h>
#include <itkLabelMap.h>
#include <itkShapeLabelObject.h>
#include <itkStatisticsLabelObject.h>
#include <itkBinaryImageToLabelMapFilter.h>
#include <itkLabelMapToLabelImageFilter.h>
#include <itkLabelImageToLabelMapFilter.h>
#include <itkShapeLabelMapFilter.h>
#include <itkStatisticsLabelMapFilter.h>
#include <itkBinaryImageToShapeLabelMapFilter.h>
#include <itkBinaryImageToStatisticsLabelMapFilter.h>

// For island removing
#include <itkBinaryShapeKeepNObjectsImageFilter.h>

#include <itkBinaryFillholeImageFilter.h>
#include <itkSliceBySliceImageFilter.h>
#include <itkSubtractImageFilter.h>
#include <itkAndImageFilter.h>

#include "itkFlatStructuringElement.h"
#include "itkBinaryErodeImageFilter.h"
#include "itkBinaryDilateImageFilter.h"

// For file operation
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>

// For threshold
#include <itkThresholdImageFilter.h>
#include <itkBinaryThresholdImageFilter.h>

#include "ITKUtils.h"

using namespace std;

int main( int argc, char *argv[] )
{
    if ( argc < 3)
    {
        cerr << "Missing Parameters " << endl;
        cerr << "Usage = " << argv[0];
        cerr << " Output InputLabelImage1 InputLabelImage2 InputLabelImage3 ..." << endl;
        return EXIT_FAILURE;
    }


    string outputName = argv[1];
    cout << "Output Name = " << outputName << endl;

    ///////// To Read the Mask images ////////////////////
    int numInputs = argc - 2;
    vector<string> inputLabelImageNames;
    vector<LabelImageType::Pointer> inputLabelImages;
    cout << "The number of input imeages = " << numInputs << endl;
    for (int i = 0; i < numInputs; i++)
    {
        inputLabelImageNames.push_back(argv[i + 2]);
        cout << "Input Label Image " << i << " Name = " << inputLabelImageNames[i] << endl;
        inputLabelImages.push_back(ReadImageFile<LabelImageType>(inputLabelImageNames[i]));
        cout << "Get the mask from Input " << i << endl;
    }

    // To check coordinate of the input label image
    SpacingType inputImageSpacing = inputLabelImages[0]->GetSpacing();
    OriginType  inputImageOrigin  = inputLabelImages[0]->GetOrigin();
    RegionType  inputImageRegion  = inputLabelImages[0]->GetLargestPossibleRegion();
    SizeType    inputImageSize    = inputImageRegion.GetSize();


    /////////////////////////////

    cout << "Input Image Spacing = " << inputImageSpacing << endl;
    cout << "Input Image Origin = " << inputImageOrigin << endl;
    cout << "Input Image Size = " << inputImageSize << endl << endl << endl;


    // Image ROI
    MaskImageType::RegionType roiRegion;
    MaskImageType::IndexType roiStart;
    MaskImageType::IndexType roiEnd;

    roiStart.Fill(1000);
    roiEnd.Fill(0);

    for (int i = 0; i < numInputs; i++)
    {
        MaskImageType::RegionType roiRegion1;
        MaskImageType::IndexType roiStart1;
        MaskImageType::IndexType roiEnd1;
    
        roiRegion1 = GetRoi(inputLabelImages[i]);
        ExpandRoi(inputLabelImages[i], roiRegion1, 10);
        RegionToIndex(roiRegion1, roiStart1, roiEnd1);

        for (unsigned i = 0; i < Dimension; i++)
        {
            if (roiStart[i] > roiStart1[i]) roiStart[i] = roiStart1[i];
            if (roiEnd[i] < roiEnd1[i]) roiEnd[i] = roiEnd1[i];
        }
    }

    roiRegion = RoiIndexToRegion(roiStart, roiEnd);

    cout << endl << "ROI start " << roiStart << endl;
    cout << "ROI end " << roiEnd << endl;
    cout << "ROI size " << roiRegion.GetSize() << endl;

    WriteROI(roiRegion,outputName);

    return EXIT_SUCCESS;
}


<|start_filename|>Tools/ContourTools/GTVs2ITV.cxx<|end_filename|>
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkResampleImageFilter.h>
#include <itkCastImageFilter.h>
#include <itkCommand.h>

#include <itkRegionOfInterestImageFilter.h>
#include <itkExtractImageFilter.h>

// Operations on Images

#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIterator.h>
#include <vnl/vnl_math.h>

// For Feature
#include <itkRescaleIntensityImageFilter.h>


#include <itkMinimumMaximumImageCalculator.h>
#include <itkStatisticsImageFilter.h>

#include <itkScalarImageToTextureFeaturesFilter.h>
#include <itkScalarImageToRunLengthFeaturesFilter.h>

#include <itkLabelGeometryImageFilter.h>
#include <itkConnectedComponentImageFilter.h>
#include <itkRelabelComponentImageFilter.h>
#include <itkLabelStatisticsImageFilter.h>

// For Label Object Representation
#include <itkLabelObject.h>
#include <itkLabelMap.h>
#include <itkShapeLabelObject.h>
#include <itkStatisticsLabelObject.h>
#include <itkBinaryImageToLabelMapFilter.h>
#include <itkLabelMapToLabelImageFilter.h>
#include <itkLabelImageToLabelMapFilter.h>
#include <itkShapeLabelMapFilter.h>
#include <itkStatisticsLabelMapFilter.h>
#include <itkBinaryImageToShapeLabelMapFilter.h>
#include <itkBinaryImageToStatisticsLabelMapFilter.h>

// For island removing
#include <itkBinaryShapeKeepNObjectsImageFilter.h>

#include <itkBinaryFillholeImageFilter.h>
#include <itkSliceBySliceImageFilter.h>
#include <itkSubtractImageFilter.h>
#include <itkAndImageFilter.h>

#include "itkFlatStructuringElement.h"
#include "itkBinaryErodeImageFilter.h"
#include "itkBinaryDilateImageFilter.h"

// For file operation
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>

// For threshold
#include <itkThresholdImageFilter.h>
#include <itkBinaryThresholdImageFilter.h>

#include "ITKUtils.h"

using namespace std;

int main( int argc, char *argv[] )
{
    if ( argc < 3)
    {
        cerr << "Missing Parameters " << endl;
        cerr << "Usage = " << argv[0];
        cerr << " Output InputLabelImage1 InputLabelImage2 InputLabelImage3 ..." << endl;
        return EXIT_FAILURE;
    }


    string outputImageName = argv[1];
    cout << "Output Name = " << outputImageName << endl;

    ///////// To Read the Mask images ////////////////////
    int numInputs = argc - 2;
    vector<string> inputLabelImageNames;
    vector<LabelImageType::Pointer> inputLabelImages;
    cout << "The number of input imeages = " << numInputs << endl;
    for (int i = 0; i < numInputs; i++)
    {
        inputLabelImageNames.push_back(argv[i + 2]);
        cout << "Input Label Image " << i << " Name = " << inputLabelImageNames[i] << endl;
        inputLabelImages.push_back(ReadImageFile<LabelImageType>(inputLabelImageNames[i]));
        cout << "Get the mask from Input " << i << endl;
    }

    // To check coordinate of the input label image
    SpacingType inputImageSpacing = inputLabelImages[0]->GetSpacing();
    OriginType  inputImageOrigin  = inputLabelImages[0]->GetOrigin();
    RegionType  inputImageRegion  = inputLabelImages[0]->GetLargestPossibleRegion();
    SizeType    inputImageSize    = inputImageRegion.GetSize();


    /////////////////////////////

    cout << "Input Image Spacing = " << inputImageSpacing << endl;
    cout << "Input Image Origin = " << inputImageOrigin << endl;
    cout << "Input Image Size = " << inputImageSize << endl << endl << endl;


    LabelImageType::Pointer outputImage = LabelImageType::New();
    outputImage->SetRegions( inputImageRegion );
    outputImage->CopyInformation( inputLabelImages[0] );
    outputImage->Allocate();
    outputImage->FillBuffer(0);

    for (int i = 0; i < numInputs; i++)
    {
        MaskImageType::RegionType roiRegion1;
        MaskImageType::IndexType roiStart1;
        MaskImageType::IndexType roiEnd1;
    
        itk::ImageRegionIterator< LabelImageType > label(inputLabelImages[i], inputLabelImages[i]->GetBufferedRegion() );
        itk::ImageRegionIterator< LabelImageType > out(outputImage, outputImage->GetBufferedRegion() );

        for (label.GoToBegin(), out.GoToBegin(); !label.IsAtEnd(); ++label, ++out)
        {
            LabelImageType::PixelType color = label.Get();
            if (color == 307)
            {
                out.Set(307);
            }
        }
    }

    ////////To write Output images /////////////////////
    try
    {
        WriteImageFile<LabelImageType>(outputImage,outputImageName);
    }
    catch (itk::ExceptionObject &excp)
    {
        cerr << "Exception thrown while writing the series " << endl;
        cerr << excp << endl;
        return EXIT_FAILURE;
    }


    return EXIT_SUCCESS;
}


<|start_filename|>Tools/ITKUtils/ITKUtils.h<|end_filename|>
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#include <itkImage.h>
#include <itkCommand.h>

#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkRegionOfInterestImageFilter.h>

#include <itkImageRegionIterator.h>
#include <itkImageRegionIteratorWithIndex.h>

#include <itkBinaryBallStructuringElement.h>
#include <itkBinaryErodeImageFilter.h>
#include <itkBinaryDilateImageFilter.h>


// std
#include <iostream>
#include <fstream>
#include <limits>


#include <csignal>

#include <string>
using namespace std;

// Definition of Data type
using InputPixelType = float;
using InternalPixelType = float;
using OutputPixelType = float;

using LabelPixelType = short;
using MaskPixelType = short;

constexpr unsigned int Dimension = 3;
constexpr unsigned int OutputDimension = 2;
constexpr MaskPixelType maskValue = 1;

using InputImageType = itk::Image < InputPixelType, Dimension >;
using InternalImageType = itk::Image < InternalPixelType, Dimension >;
using OutputImageType = itk::Image < OutputPixelType, OutputDimension >;
using InputImage2DType = itk::Image < InputPixelType, Dimension - 1 >;

using LabelImageType = itk::Image < LabelPixelType, Dimension>;
using MaskImageType = itk::Image < MaskPixelType, Dimension>;
using MaskImage2DType = itk::Image < MaskPixelType, Dimension - 1 >;


// Definition of Properties
using SpacingType = InputImageType::SpacingType;
using OriginType = InputImageType::PointType;
using RegionType = InputImageType::RegionType;
using SizeType = InputImageType::SizeType;
using IndexType = InputImageType::IndexType;
using DirectionType = InputImageType::DirectionType;


//Definition of Morphological operation sturcture elements
using StructuringElementType = itk::BinaryBallStructuringElement< MaskPixelType, Dimension >;
using DilateFilterType = itk::BinaryDilateImageFilter<MaskImageType, MaskImageType, StructuringElementType>;
using ErodeFilterType = itk::BinaryErodeImageFilter<MaskImageType, MaskImageType, StructuringElementType>;


template <typename TImageType>
typename TImageType::Pointer ReadImageFile(string fileName)
{
  using ImageReaderType = itk::ImageFileReader<TImageType>;

  typename ImageReaderType::Pointer ImageReader = ImageReaderType::New();
  ImageReader->SetFileName(fileName);
  ImageReader->Update();

  itk::SmartPointer<TImageType> Image = ImageReader->GetOutput();

  return Image;
}


template <typename TImageType>
void WriteImageFile(typename TImageType::Pointer outputImage, string fileName)
{
  ////////To write Output images /////////////////////
  using ImageWriterType = itk::ImageFileWriter<TImageType>;

  typename ImageWriterType::Pointer OutImageWriter = ImageWriterType::New();
  OutImageWriter->SetUseCompression(true);
  OutImageWriter->SetInput(outputImage);
  OutImageWriter->SetFileName( fileName );

  try
  {
    OutImageWriter->Update();
  }
  catch (itk::ExceptionObject &excp)
  {
    std::cerr << "Exception thrown while writing the series " << std::endl;
    std::cerr << excp << std::endl;
    //return EXIT_FAILURE;
  }
}


template <typename TImageType>
typename TImageType::Pointer ApplyRoi(typename TImageType::Pointer inputImage, itk::ImageRegion<Dimension> inputRegion)
{
  typename TImageType::Pointer roiImage;

  using RoiFilterType = itk::RegionOfInterestImageFilter<TImageType, TImageType>;
  typename RoiFilterType::Pointer RoiFilter = RoiFilterType::New();
  RoiFilter->SetRegionOfInterest(inputRegion);

  RoiFilter->SetInput( inputImage );
  RoiFilter->Update();

  roiImage = RoiFilter->GetOutput();

  return roiImage;
}

MaskImageType::Pointer unionImages(vector<MaskImageType::Pointer> inputMaskImages);

namespace FGC
{

/*
* createImage
*/

template<typename itkImageType>
typename itkImageType::Pointer createImage(typename itkImageType::SizeType size, int iniValue)
{

  typename itkImageType::Pointer img = itkImageType::New();
  typename itkImageType::SpacingType spacing;
  typename itkImageType::IndexType start;
  start.Fill(0);
  spacing.Fill(1);
  spacing[itkImageType::ImageDimension - 1] = 2;

  typename itkImageType::RegionType region;
  img->SetSpacing(spacing);
  region.SetSize(size);
  region.SetIndex(start);
  img->SetRegions(region);
  img->Allocate();
  img->FillBuffer(static_cast<typename itkImageType::PixelType>(iniValue));

  const typename itkImageType::SpacingType &sp = img->GetSpacing();
  std::cout << "Spacing = ";
  std::cout << sp[0] << ", " << sp[1] << std::endl;


  return img;

}


template<typename itkImageType>
void FindITKImageROI(typename itkImageType::Pointer &im, std::vector<long> &imROI)
{

  typename itkImageType::IndexType roiStart;
  typename itkImageType::IndexType roiEnd;
  typename itkImageType::IndexType start;
  typename itkImageType::SizeType size;

  size = im->GetLargestPossibleRegion().GetSize();
  start = im->GetLargestPossibleRegion().GetIndex();


  roiStart[0] = 0; roiStart[1] = 0; roiStart[2] = 0;
  roiEnd[0] = 0; roiEnd[1] = 0; roiEnd[2] = 0;

  unsigned int ndims = im->GetImageDimension();

  bool foundLabel = false;
  itk::ImageRegionIteratorWithIndex< itkImageType > label(im, im->GetBufferedRegion() );
  for (label.GoToBegin(); !label.IsAtEnd(); ++label)
  {
    if (label.Get() != 0)
    {
      typename itkImageType::IndexType idx = label.GetIndex();
      for (unsigned i = 0; i < ndims; i++)
      {
        if (!foundLabel)
        {
          roiStart[i] = idx[i];
          roiEnd[i] = idx[i];
        }
        else
        {
          if (idx[i] <= roiStart[i])
          {
            roiStart[i] = idx[i];
          }
          if (idx[i] >= roiEnd[i])
          {
            roiEnd[i] = idx[i];
          }
        }
      }
      foundLabel = true;
    }
  }

  int radius = 17;
  for (unsigned i = 0; i < ndims; i++)
  {
    int diff = static_cast< int > (roiStart[i] - radius);
    if (diff >= start[i])
    {
      roiStart[i] -= radius;
    }
    else
    {
      roiStart[i] = start[i];
    }
    roiEnd[i] = (static_cast<unsigned int>(roiEnd[i] + radius) < size[i]) ? (roiEnd[i] + radius) : size[i] - 1;

  }

  // copy ROI to vector
  imROI.resize(6);
  for (unsigned i = 0; i < 3; i++)
  {
    imROI[i] = roiStart[i];
    imROI[i + 3] = roiEnd[i];
  }
}

template<typename PixelType>
void ExtractITKImageROI(typename itk::Image<PixelType, 3>::Pointer  im, const std::vector<long> &imROI, \
                        std::vector<PixelType> &imROIVec)
{

  // Copy itk image ROI to vector
  using ImageType = itk::Image<PixelType, 3>;
  typename ImageType::IndexType index;
  long i, j, k, kk, DIMXYZ;

  DIMXYZ = (imROI[3] - imROI[0]) * (imROI[4] - imROI[1]) * (imROI[5] - imROI[2]);
  imROIVec.clear();
  imROIVec.resize(DIMXYZ);
  kk = 0;
  for (k = imROI[2]; k < imROI[5]; k++)
    for (j = imROI[1]; j < imROI[4]; j++)
      for (i = imROI[0]; i < imROI[3]; i++)
      {
        index[0] = i; index[1] = j; index[2] = k;
        imROIVec[kk++] = im->GetPixel(index);
      }
}

template<typename PixelType>
void UpdateITKImageROI(const std::vector<PixelType> &imROIVec, const std::vector<long> &imROI,  \
                       typename itk::Image<PixelType, 3>::Pointer im)
{

  using ImageType = itk::Image<PixelType, 3>;
  typename ImageType::IndexType index;
  long i, j, k, kk;

  // Set non-ROI as zeros
  im->FillBuffer(0);
  kk = 0;
  for (k = imROI[2]; k < imROI[5]; k++)
    for (j = imROI[1]; j < imROI[4]; j++)
      for (i = imROI[0]; i < imROI[3]; i++)
      {
        index[0] = i; index[1] = j; index[2] = k;
        im->SetPixel(index, imROIVec[kk++]);
      }
}




/**
 * readImage
 */
template< typename TImageType >
typename TImageType::Pointer readImage(const char *fileName)
{
  using ImageReaderType = itk::ImageFileReader< TImageType >;
  typename ImageReaderType::Pointer reader = ImageReaderType::New();
  reader->SetFileName(fileName);

  typename TImageType::Pointer image;

  try
  {
    reader->Update();
    image = reader->GetOutput();
  }
  catch ( itk::ExceptionObject &err)
  {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    raise(SIGABRT);
  }

  return image;
}

/**
 * writeImage
 */
template< typename TImageType > void writeImage(typename TImageType::Pointer img, const char *fileName)
{
  using WriterType = itk::ImageFileWriter< TImageType >;

  typename WriterType::Pointer writer = WriterType::New();
  writer->SetFileName( fileName );
  writer->SetInput(img);
  writer->UseCompressionOn();

  try
  {
    writer->Update();
  }
  catch ( itk::ExceptionObject &err)
  {
    std::cout << "ExceptionObject caught !" << std::endl;
    std::cout << err << std::endl;
    raise(SIGABRT);
  }
}

template<typename VType>
void WriteVectorIntoFile(const char *fnSave, const std::vector<VType> &vec)
{

  try
  {
    std::ofstream outfile(fnSave);
    outfile.write((const char *)&vec[0], vec.size()*sizeof(VType));
  }
  catch ( itk::ExceptionObject &err)
  {
    std::cout << "Fail to write to file " << fnSave << std::endl;
    std::cout << err << std::endl;
    raise(SIGABRT);
  }
}

template<typename VType>
void LoadVectorIntoFile(const char *fnLoad, std::vector<VType> &vec, const long VEC_SIZE)
{

  try
  {
    std::ifstream infile(fnLoad);
    vec.resize(VEC_SIZE);
    infile.read((char *)&vec[0], vec.size()*sizeof(VType));
  }
  catch ( itk::ExceptionObject &err)
  {
    std::cout << "Fail to load file " << fnLoad << std::endl;
    std::cout << err << std::endl;
    raise(SIGABRT);
  }
}
}// douher

MaskImageType::Pointer GetMaskImage(LabelImageType::Pointer);
MaskImageType::Pointer GetMaskImage(LabelImageType::Pointer, LabelPixelType);

MaskImageType::RegionType RoiIndexToRegion(MaskImageType::IndexType, MaskImageType::IndexType);
void RegionToIndex(MaskImageType::RegionType, MaskImageType::IndexType &, MaskImageType::IndexType &);

MaskImageType::RegionType GetRoi(MaskImageType::Pointer);
void ExpandRoi(MaskImageType::Pointer, MaskImageType::RegionType &, const double objectSize = 10);
void BoundingCheck(MaskImageType::Pointer, InputImageType::Pointer, MaskImageType::RegionType &, InputImageType::RegionType &);
void BoundingCheck(MaskImageType::Pointer, LabelImageType::Pointer, MaskImageType::RegionType &, InputImageType::RegionType &);

MaskImageType::RegionType ReadROI(string);
void WriteROI(MaskImageType::RegionType, string);

// To print the progress
class ShowProgressObject
{
public:
  ShowProgressObject( itk::ProcessObject *o )
  {
    m_Process = o;
  }
  void ShowProgress()
  {
    std::cout << "\rProgress: "
              << static_cast<unsigned int>( 100.0 * m_Process->GetProgress() ) << "%";
  }
  itk::ProcessObject::Pointer m_Process;
};


<|start_filename|>Tools/ITKUtils/itkFastGrowCutSegmentationImageFilter.h<|end_filename|>
#ifndef __itkFastGrowCutSegmentationImageFilter_h
#define __itkFastGrowCutSegmentationImageFilter_h

#include "itkImage.h"
#include "itkImageToImageFilter.h"
#include "itkSimpleDataObjectDecorator.h"
#include "itkVectorContainer.h"
//#include "itkCommand.h"

#include "vcl_list.h"

#include <vcl_vector.h>

#include "FastGrowCutSegmenter.h"

//#include "itkGrowCutSegmentationUpdateFilter.h"

namespace itk
{

/**  \class GrowCutSegmentationFilter
 * \brief Segmentation or multiple objects based on a set of gestures.
 *
 * Given a mask image containing the gestures for foreground classes
 * and their background, employs grow cut segmentation to produce
 * segmentation. Supports passing an existing segmentation with
 * additional gestures for editing a segmentation produced as a
 * result of the same or different algorithm.
 *
 * The filter is based on the paper "GrowCut:Interactive Multi-Label
 * N-D Image Segmentation By Cellular Automata", <NAME>,
 * <NAME>
 *
 * Modified Version: The inputs consist of the input intensity image,
 * optionally the input Label Image, and the Seed image.  The Label
 * image and the Seed image default to using unsigned chars and
 * double respectively.
 *
 * This algorithm is implemented scalar images. Vector Images are not
 * supported.
 *
 *
**/


/* template<class TInputImage,  */
/*   class TOutputImage, class TLabelPixelType = short,  */
/*   class TSeedPixelType = float >  */
template<class TInputImage,
         class TOutputImage>
class FastGrowCutSegmentationImageFilter: public ImageToImageFilter<TInputImage, TOutputImage>
{

public:
    /** Standard class type alias. */
    using Self = FastGrowCutSegmentationImageFilter;
    using Superclass = ImageToImageFilter<TInputImage, TOutputImage>;
    using Pointer = SmartPointer<Self>;
    using ConstPointer = SmartPointer<const Self>;

    /** Method for creation through the object factory. */
    itkNewMacro(Self);

    /** Run-time type information (and related methods).  */
    itkTypeMacro(FastGrowCutSegmentationImageFilter,
                 ImageToImageFilter);

    /** Image related type alias. */
    itkStaticConstMacro(ImageDimension, unsigned int,
                        TInputImage::ImageDimension );

    using InputImageType = TInputImage;
    using InputImagePointer = typename InputImageType::Pointer;
    using InputImageConstPointer = typename InputImageType::ConstPointer;

    using InputPixelType = typename InputImageType::PixelType;
    using InputIndexType = typename InputImageType::IndexType;
    using SizeType = typename InputImageType::SizeType;

    using OutputImageType = TOutputImage;
    using OutputImagePointer = typename OutputImageType::Pointer;
    using OutputImageRegionType = typename OutputImageType::RegionType;
    using OutputPixelType = typename OutputImageType::PixelType;
    using OutputIndexType = typename OutputImageType::IndexType;
    using OutputSizeType = typename InputImageType::SizeType;


    /** Smart Pointer type to a DataObject. */
    using DataObjectPointer = typename DataObject::Pointer;


    /** Image dimension constants */
    itkStaticConstMacro(InputImageDimension, unsigned int,
                        TInputImage::ImageDimension);
    itkStaticConstMacro(OutputImageDimension, unsigned int,
                        TOutputImage::ImageDimension);

    /** Index type alias support. */
    using IndexType = Index<itkGetStaticConstMacro(InputImageDimension)>;

    /** InputSizeType type alias support **/
    using InputSizeType = typename InputImageType::SizeType;

    /** SeedImage type alias support */
    /* indicates the Seed of a label for a given cell */
    using SeedImageType = TOutputImage;

    /** SeedImagePointer type alias support. */
    using SeedImagePointer = typename SeedImageType::Pointer;

    /** Set the Input Image **/
    void SetInputImage( const InputImageType *in)
    {
        this->ProcessObject::SetNthInput(0, const_cast< InputImageType *>(in) );
    }

    const InputImagePointer GetInputImage( );

    /** Set/Get the Label Image **/
    void SetLabelImage( const OutputImageType *f)
    {
        this->ProcessObject::SetNthInput(2, const_cast< OutputImageType *>(f) );
        m_LabelImage = static_cast< OutputImageType *>(this->ProcessObject::GetInput(2));
    }

    const OutputImagePointer GetLabelImage();

    /** Get the Seed image **/
    const SeedImagePointer GetSeedImage();

    /** Set the Seed Image **/
    void SetSeedImage( const SeedImageType *s )
    {
        this->ProcessObject::SetNthInput(1, const_cast< SeedImageType *>(s));
        //m_SeedImage = static_cast< SeedImageType *>(this->ProcessObject::GetInput(2));
    }

    const SeedImagePointer GetUpdatedSeedImage();

    itkSetMacro(InitializationFlag, bool);

    //processing functions
    void Initialization();
    void RunFGC();

protected:

    FastGrowCutSegmentationImageFilter();
    ~FastGrowCutSegmentationImageFilter() {};

    // Override since the filter needs all the data for the algorithm
    void GenerateInputRequestedRegion();

    // Override since the filter produces the entire dataset
    void EnlargeOutputRequestedRegion(DataObject *output);

    void GenerateData();

    void Initialize(OutputImageType *output);

    void PrintSelf ( std::ostream &os, Indent indent ) const;

    void GrowCutSlowROI( TOutputImage *);


private:

    FastGrowCutSegmentationImageFilter(const Self &); //purposely not implemented
    void operator=(const Self &); // purposely not implemented

    void GetRegionOfInterest();

    void ComputeLabelVolumes(TOutputImage *outputImage, vcl_vector< unsigned > &volumes, vcl_vector< unsigned > &phyVolumes);

    OutputImagePointer                       m_LabelImage;
    SeedImagePointer                         m_SeedImage;


    std::vector<OutputPixelType> m_imSeedVec;
    std::vector<OutputPixelType> m_imLabVec;
    std::vector<InputPixelType> m_imSrcVec;
    std::vector<long> m_imROI;

    //logic code
    FGC::FastGrowCut<InputPixelType, OutputPixelType> *m_fastGC;

    //state variables
    bool m_InitializationFlag;

};

} // namespace itk


#ifndef ITK_MANUAL_INSTANTIATION
#include "itkFastGrowCutSegmentationImageFilter.hxx"
#endif

#endif


<|start_filename|>Tools/ContourTools/ROICropImage.cxx<|end_filename|>
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkResampleImageFilter.h>
#include <itkCastImageFilter.h>
#include <itkCommand.h>

#include <itkRegionOfInterestImageFilter.h>
#include <itkExtractImageFilter.h>

// Operations on Images

#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIterator.h>
#include <vnl/vnl_math.h>

// For Feature
#include <itkRescaleIntensityImageFilter.h>


#include <itkMinimumMaximumImageCalculator.h>
#include <itkStatisticsImageFilter.h>

#include <itkScalarImageToTextureFeaturesFilter.h>
#include <itkScalarImageToRunLengthFeaturesFilter.h>

#include <itkLabelGeometryImageFilter.h>
#include <itkConnectedComponentImageFilter.h>
#include <itkRelabelComponentImageFilter.h>
#include <itkLabelStatisticsImageFilter.h>

// For Label Object Representation
#include <itkLabelObject.h>
#include <itkLabelMap.h>
#include <itkShapeLabelObject.h>
#include <itkStatisticsLabelObject.h>
#include <itkBinaryImageToLabelMapFilter.h>
#include <itkLabelMapToLabelImageFilter.h>
#include <itkLabelImageToLabelMapFilter.h>
#include <itkShapeLabelMapFilter.h>
#include <itkStatisticsLabelMapFilter.h>
#include <itkBinaryImageToShapeLabelMapFilter.h>
#include <itkBinaryImageToStatisticsLabelMapFilter.h>

// For island removing
#include <itkBinaryShapeKeepNObjectsImageFilter.h>

#include <itkBinaryFillholeImageFilter.h>
#include <itkSliceBySliceImageFilter.h>
#include <itkSubtractImageFilter.h>
#include <itkAndImageFilter.h>

#include "itkFlatStructuringElement.h"
#include "itkBinaryErodeImageFilter.h"
#include "itkBinaryDilateImageFilter.h"

// For file operation
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>

// For threshold
#include <itkThresholdImageFilter.h>
#include <itkBinaryThresholdImageFilter.h>

#include "ITKUtils.h"

using namespace std;

int main( int argc, char *argv[] )
{
    if ( argc < 3)
    {
        cerr << "Missing Parameters " << endl;
        cerr << "Usage = " << argv[0];
        cerr << " OutputLabelImage inputImage [ROIfile or start(x,y,z) size(x,y,z)]" << endl;
        return EXIT_FAILURE;
    }


    string outputImageName = argv[1];
    string inputImageName = argv[2];
    string roiName = argv[3];
    
    cout << "Output Label Image Name = " << outputImageName << endl;
    cout << "Input Image 1 Name = " << inputImageName << endl;
    
    ///////// To Read the Reference image ////////////////////
    cout << "Get the Input Image" << endl;
    InputImageType::Pointer inputImage = ReadImageFile<InputImageType>(inputImageName);
    

    // To check coordinate of the input label image
    SpacingType inputImageSpacing = inputImage->GetSpacing();
    OriginType  inputImageOrigin  = inputImage->GetOrigin();
    RegionType  inputImageRegion  = inputImage->GetLargestPossibleRegion();
    SizeType    inputImageSize    = inputImageRegion.GetSize();


    /////////////////////////////

    cout << "Input Image Spacing = " << inputImageSpacing << endl;
    cout << "Input Image Origin = " << inputImageOrigin << endl;
    cout << "Input Image Size = " << inputImageSize << endl << endl << endl;


    // Read ROI
    InputImageType::RegionType roiRegion;
    InputImageType::IndexType roiStart;
    InputImageType::IndexType roiEnd;
    InputImageType::Pointer outputImage;

    roiRegion = ReadROI(roiName);
    cout << roiRegion << endl;
    RegionToIndex(roiRegion, roiStart, roiEnd);
    outputImage = ApplyRoi<InputImageType>(inputImage, roiRegion);

    cout << "writeImage_spacing = " << outputImage->GetSpacing() << endl ;
    cout << "writeImage_origin  = "  << outputImage->GetOrigin() << endl ;
    cout << "writeImage_LargestPossibleRegion = " << outputImage->GetLargestPossibleRegion() << endl ;


    ////////To write Output images /////////////////////
    WriteImageFile<InputImageType>(outputImage,outputImageName);


    std::cout << "Saved the normal tissue and tumor boundary image" << std::endl;

    return EXIT_SUCCESS;
}


<|start_filename|>Tools/DICOMTools/DICOMTagReader.cxx<|end_filename|>
#include "gdcmReader.h"
#include "gdcmAttribute.h"

#include "gdcmGlobal.h"
#include "gdcmDicts.h"
#include "gdcmDict.h"
#include "gdcmCSAHeader.h"
#include "gdcmPrivateTag.h"



void recursive_read(const gdcm::DataSet& ds)
{
    static int depth = 0;
    const gdcm::Global& g = gdcm::Global::GetInstance(); // sum of all knowledge !
    const gdcm::Dicts &dicts = g.GetDicts();
    const gdcm::Dict &dict = dicts.GetPublicDict();

    for (gdcm::DataSet::ConstIterator it = ds.Begin(); it!=ds.End(); ++it)
    {
        const gdcm::DataElement& elem = *it;

        if(!elem.IsEmpty())
        {    
            const gdcm::Tag& tag = elem.GetTag();

            for(int d = 0; d < depth; d++)
                std::cout << ' ';
            std::cout << dict.GetDictEntry(tag).GetKeyword() << ": ";
            std::cout << elem.GetValue() << std::endl;
			try
			{
				gdcm::SmartPointer<gdcm::SequenceOfItems> ssqi = elem.GetValueAsSQ();
				if (ssqi != NULL)
				{
					for (int i = 1; i <= ssqi->GetNumberOfItems(); i++)
					{
						const gdcm::Item & item = ssqi->GetItem(i); // Item start at #1
						const gdcm::DataSet& nested_ds = item.GetNestedDataSet();

						depth++;
						recursive_read(nested_ds);
						depth--;
					}
				}
            }
            catch (const gdcm::Exception& e)
            {
                std::cerr << e.what() << std::endl;
            }
        }
     }
}

int main( int argc, char* argv[] )
{
    if( argc < 2 )
    {
        std::cerr << "Usage: " << argv[0] << " DicomFilename " << std::endl;
        return EXIT_FAILURE;
    }

    gdcm::Reader reader;
    reader.SetFileName(argv[1]);
    if( !reader.Read() )
    {
        std::cerr << "Cannot read DICOM file!!" << std::endl;
        return 1;
    }

    gdcm::MediaStorage ms;
    ms.SetFromFile( reader.GetFile() );
    if( ms != gdcm::MediaStorage::SpacialRegistrationStorage)
    {
        std::cerr << "No DICOM REG!!" << std::endl;
    }

    const gdcm::DataSet& ds = reader.GetFile().GetDataSet();
    recursive_read(ds);
    

    return EXIT_SUCCESS;
}


<|start_filename|>Tools/FeatureExtraction/RunLengthFeatures.cxx<|end_filename|>
#include <iostream>
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include <itkDenseFrequencyContainer2.h>
#include "itkRegionOfInterestImageFilter.h"
#include "itkScalarImageToRunLengthFeaturesFilter.h"
#include "itkScalarImageToTextureFeaturesFilter.h"
#include "itkOffset.h"
#include "itkImageBase.h"
#include "itkIndex.h"
#include <cmath>

//definitions of used types
using InternalImageType = itk::Image<float, 3>;
using NeighborhoodType = itk::Neighborhood<float, 3>;
using Image2runLengthType = itk::Statistics::ScalarImageToRunLengthMatrixFilter<InternalImageType>;
using HistogramType = Image2runLengthType::HistogramType;
using Hist2FeaturesType = itk::Statistics::HistogramToRunLengthFeaturesFilter<HistogramType>;
using HistogramType = Image2runLengthType::HistogramType;
using OffsetType = InternalImageType::OffsetType;
using SpacingType = InternalImageType::SpacingType;

//calculate features for one offset
float * calcTextureFeatureImage (OffsetType offset,
    InternalImageType::Pointer inputImage, int numBins,
	float minVal, float maxVal, InternalImageType::Pointer maskImage)

{
    
	// Initialize an array of floats to store global features
	static float featuresV[10];

	// create an rlm generator
    Image2runLengthType::Pointer rlGenerator = Image2runLengthType::New();

	// define offset (hard-coded for now). should assign from the input offset.
	OffsetType list_offset;
	list_offset[0] = 1;
	list_offset[1] = 0;
	list_offset[2] = 0;

	// set offset
	rlGenerator->SetOffset(offset);
    rlGenerator->SetNumberOfBinsPerAxis(numBins); //reasonable number of bins
    rlGenerator->SetPixelValueMinMax(minVal, maxVal); //for input UCHAR pixel type
    
	/*
	InternalImageType::PointType pointMin;
    InternalImageType::PointType pointMax;
	itk::Index<3> minIndex;
	itk::Index<3> maxIndex;
	minIndex.SetElement(0,0);
	minIndex.SetElement(1,0);
	minIndex.SetElement(2,0);
	maxIndex.SetElement(0,64);
	maxIndex.SetElement(1,64);
	maxIndex.SetElement(2,64);	

    maskImage->TransformIndexToPhysicalPoint(minIndex, pointMin);
    maskImage->TransformIndexToPhysicalPoint(maxIndex, pointMax);

    double distMax = pointMin.EuclideanDistanceTo(pointMax);

	std::cout << "----- dist ------ ";
	std::cout << distMax << std::endl;

	rlGenerator -> SetDistanceValueMinMax (0,distMax);
	*/

	// create pointer to histogram object
	Hist2FeaturesType::Pointer featureCalc = Hist2FeaturesType::New();	

	// assign image and mask to the rlm generator
	rlGenerator -> SetInput(inputImage);
	rlGenerator -> SetMaskImage(maskImage);
	
	featureCalc -> SetInput (rlGenerator->GetOutput());
    featureCalc -> Update();

	// set correct max distance for this direction
	double numLevels = rlGenerator -> GetMax();	
	double voxelLength = 0;
	SpacingType spacing = inputImage -> GetSpacing();
	if (abs(offset[0]) == 1){
		voxelLength = spacing[0];
	}
	if (abs(offset[1]) == 1){
		voxelLength = pow(voxelLength,2) + pow(spacing[1],2);
		voxelLength = sqrt(voxelLength);
	}
	if (abs(offset[2]) == 1){
		voxelLength = pow(voxelLength,2) + pow(spacing[2],2);
		voxelLength = sqrt(voxelLength);
	}
	double dMax = voxelLength * numLevels;
	std::cout << "Maximum distance for: ";
	std::cout << offset << std::endl;
	std::cout << " is: ";
	std::cout << dMax << std::endl;
	rlGenerator -> SetDistanceValueMinMax (0,dMax);
	featureCalc -> Update();

	featuresV[0] = featureCalc->GetFeature(Hist2FeaturesType::ShortRunEmphasis); 
	featuresV[1] = featureCalc->GetFeature(Hist2FeaturesType::LongRunEmphasis); 
	featuresV[2] = featureCalc->GetFeature(Hist2FeaturesType::GreyLevelNonuniformity); 
	featuresV[3] = featureCalc->GetFeature(Hist2FeaturesType::RunLengthNonuniformity); 
	featuresV[4] = featureCalc->GetFeature(Hist2FeaturesType::LowGreyLevelRunEmphasis); 
	featuresV[5] = featureCalc->GetFeature(Hist2FeaturesType::HighGreyLevelRunEmphasis); 
	featuresV[6] = featureCalc->GetFeature(Hist2FeaturesType::ShortRunLowGreyLevelEmphasis); 
	featuresV[7] = featureCalc->GetFeature(Hist2FeaturesType::ShortRunHighGreyLevelEmphasis); 
	featuresV[8] = featureCalc->GetFeature(Hist2FeaturesType::LongRunLowGreyLevelEmphasis); 
	featuresV[9] = featureCalc->GetFeature(Hist2FeaturesType::LongRunHighGreyLevelEmphasis); 
    
	return featuresV;

}

int main(int argc, char*argv[])
{
  if(argc < 5)
    {
    std::cerr << "Usage: " << argv[0] << " Required image.mha" << std::endl;
    return EXIT_FAILURE;
    }
  
  std::string fileName = argv[1];
  std::string numBinsStr = argv[2];
  int numBins = atoi(numBinsStr.c_str());
  std::string minValStr = argv[3];
  float minVal = atof(minValStr.c_str());
  std::string maxValStr = argv[4];
  float maxVal = atof(maxValStr.c_str());
  std::string maskFileName = argv[5];

  // read image
  using ReaderType = itk::ImageFileReader<InternalImageType>;
  ReaderType::Pointer reader=ReaderType::New();
  reader->SetFileName(fileName);
  reader->Update();
  InternalImageType::Pointer image=reader->GetOutput();

  // read mask
  ReaderType::Pointer maskReader=ReaderType::New();
  maskReader->SetFileName(maskFileName);
  maskReader->Update();
  InternalImageType::Pointer mask=maskReader->GetOutput();

  // neighborhood to get offset vectors
  NeighborhoodType neighborhood;
  neighborhood.SetRadius(1);
  unsigned int centerIndex = neighborhood.GetCenterNeighborhoodIndex();
  
  OffsetType offset;

  using WriterType = itk::ImageFileWriter<InternalImageType>;
  WriterType::Pointer writer=WriterType::New();

  for ( unsigned int d = 0; d < centerIndex; d++ )
  {
      offset = neighborhood.GetOffset(d);

	  float * featuresV;

	  featuresV = calcTextureFeatureImage(offset, image, numBins, minVal, maxVal, mask);

	  for ( int i = 0; i < 10; i++ )
	  {
		  std::cout << "*(featuresV + " << i << ") : ";
		  std::cout << *(featuresV + i) << std::endl;
      }
	
  }

  return EXIT_SUCCESS;
}


<|start_filename|>Tools/ContourTools/STAPLEComparison.cxx<|end_filename|>
#if defined(_MSC_VER)
#pragma warning(disable : 4786)
#endif

#include <itkCastImageFilter.h>
#include <itkCommand.h>
#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkResampleImageFilter.h>

#include <itkPasteImageFilter.h>
#include <itkCropImageFilter.h>
#include <itkRegionOfInterestImageFilter.h>

// Operations on Images

#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIterator.h>
#include <vnl/vnl_math.h>

// For Feature
#include <itkRescaleIntensityImageFilter.h>

#include <itkMinimumMaximumImageCalculator.h>
#include <itkStatisticsImageFilter.h>

#include <itkScalarImageToRunLengthFeaturesFilter.h>
#include <itkScalarImageToTextureFeaturesFilter.h>

#include <itkConnectedComponentImageFilter.h>
#include <itkLabelGeometryImageFilter.h>
#include <itkLabelStatisticsImageFilter.h>
#include <itkRelabelComponentImageFilter.h>

// For Label Object Representation
#include <itkBinaryImageToLabelMapFilter.h>
#include <itkBinaryImageToShapeLabelMapFilter.h>
#include <itkBinaryImageToStatisticsLabelMapFilter.h>
#include <itkLabelImageToLabelMapFilter.h>
#include <itkLabelMap.h>
#include <itkLabelMapToLabelImageFilter.h>
#include <itkLabelObject.h>
#include <itkShapeLabelMapFilter.h>
#include <itkShapeLabelObject.h>
#include <itkStatisticsLabelMapFilter.h>
#include <itkStatisticsLabelObject.h>

// For island removing
#include <itkBinaryShapeKeepNObjectsImageFilter.h>

#include <itkAndImageFilter.h>
#include <itkBinaryFillholeImageFilter.h>
#include <itkSliceBySliceImageFilter.h>
#include <itkSubtractImageFilter.h>

#include "itkBinaryDilateImageFilter.h"
#include "itkBinaryErodeImageFilter.h"
#include "itkFlatStructuringElement.h"

#include <itkMaskedSTAPLEImageFilter.h>
#include <itkSTAPLEImageFilter.h>

// For file operation
#include <fstream>
#include <iomanip>
#include <iostream>
#include <string>

// For threshold
#include <itkBinaryThresholdImageFilter.h>
#include <itkThresholdImageFilter.h>

using namespace std;

#include "ITKUtils.h"

int main(int argc, char *argv[])
{
    if (argc < 4)
    {
        cerr << "Missing Parameters " << endl;
        cerr << "Usage = " << argv[0];
        cerr << " OutputImage OutputCsv Label InputLabelImage1 InputLabelImage2 "
                "InputLabelImage3 ..."
             << endl;
        return EXIT_FAILURE;
    }

    int numInputs = argc - 4;

    cout << "The number of input imeages = " << numInputs << endl;

    string outputImageName = argv[1];
    string outputCsvFileName = argv[2];

    // Output Text file open
    ofstream outfile((outputCsvFileName).c_str(), ios::out);
    if (!outfile)
    {
        cerr << "Can't open output file " << outputCsvFileName << endl;
        exit(1);
    }

    LabelPixelType selectedLabelValue = atoi(argv[3]);

    vector<string> inputLabelImageNames;
    cout << "Output Image Name = " << outputImageName << endl;
    for (int i = 0; i < numInputs; i++)
    {
        inputLabelImageNames.push_back(argv[i + 4]);
        cout << "Input Label Image " << i << " Name = " << inputLabelImageNames[i]
             << endl;
    }

    ///////// To Read the Mask image ////////////////////
    OriginType newOrigin;
    newOrigin[0] = 10000;
    newOrigin[1] = 10000;
    newOrigin[2] = 10000;
    OriginType newExtent;
    newExtent[0] = -10000;
    newExtent[1] = -10000;
    newExtent[2] = -10000;
    SizeType newSize;
    RegionType newRegion;

    vector<MaskImageType::Pointer> inputMaskImages;
    for (int i = 0; i < numInputs; i++)
    {
        MaskImageType::Pointer maskImage;
        if (selectedLabelValue > 0)
        {
            maskImage =
                GetMaskImage(ReadImageFile<LabelImageType>(inputLabelImageNames[i]),
                             selectedLabelValue);
        }
        else
        {
            maskImage =
                GetMaskImage(ReadImageFile<LabelImageType>(inputLabelImageNames[i]));
        }

        auto origin = maskImage->GetOrigin();
        auto spacing = maskImage->GetSpacing();
        auto size = maskImage->GetLargestPossibleRegion().GetSize();
        inputMaskImages.push_back(maskImage);
        for (int j = 0; j < 3; j++)
        {
            auto extent = origin[j] + size[j] * spacing[j];
            if (origin[j] < newOrigin[j])
                newOrigin[j] = origin[j];
            if (extent > newExtent[j])
                newExtent[j] = extent;
            newSize[j] = (newExtent[j] - newOrigin[j]) / spacing[j];
        }
    }
    newRegion.SetIndex({{0, 0, 0}});
    newRegion.SetSize(newSize);
    cout << newOrigin << newExtent << newSize << endl;

    for (int i = 0; i < inputMaskImages.size(); i++)
    {
        auto oldMaskImage = inputMaskImages[i];
        auto newMaskImage = MaskImageType::New();
        newMaskImage->CopyInformation(oldMaskImage);
        newMaskImage->SetOrigin(newOrigin);
        newMaskImage->SetRegions(newRegion);
        newMaskImage->Allocate();
        newMaskImage->FillBuffer(0);

        // To check coordinate of the input label image
        auto oldMaskSpacing = oldMaskImage->GetSpacing();
        auto oldMaskOrigin = oldMaskImage->GetOrigin();
        auto oldMaskRegion = oldMaskImage->GetLargestPossibleRegion();
        auto oldMaskSize = oldMaskRegion.GetSize();

        /////////////////////////////
        cout << "Old Mask Spacing = " << oldMaskSpacing << endl;
        cout << "Old Mask Origin = " << oldMaskOrigin << endl;
        cout << "Old Mask  Size = " << oldMaskSize << endl
             << endl;

        OriginType diffOrigin;
        MaskImageType::IndexType destinationIndex;
        for (int j = 0; j < 3; j++)
        {
            diffOrigin[j] = oldMaskOrigin[j] - newOrigin[j];
            destinationIndex[j] = diffOrigin[j] / oldMaskSpacing[j];
        }
        cout << diffOrigin << destinationIndex << endl;

        using PasteFilterType = itk::PasteImageFilter<MaskImageType>;
        auto pasteFilter = PasteFilterType::New();
        pasteFilter->SetSourceImage(oldMaskImage);
        pasteFilter->SetDestinationImage(newMaskImage);
        pasteFilter->SetSourceRegion(oldMaskRegion);
        pasteFilter->SetDestinationIndex(destinationIndex);
        pasteFilter->Update();

        inputMaskImages[i] = pasteFilter->GetOutput();

        // To check coordinate of the label image
        auto newMaskOrigin = inputMaskImages[i]->GetOrigin();
        auto newMaskRegion = inputMaskImages[i]->GetBufferedRegion();
        auto newMaskSize = newMaskRegion.GetSize();

        /////////////////////////////
        cout << "New Mask Origin = " << newMaskOrigin << endl;
        cout << "New Mask Size = " << newMaskSize << endl
             << endl;
        //WriteImageFile<MaskImageType>(inputMaskImages[i], inputLabelImageNames[i]+"test-label.nrrd");
    }
    cout << "Get the mask from Input" << endl;
    MaskImageType::Pointer unionMaskImage = unionImages(inputMaskImages);
    //WriteImageFile<MaskImageType>(unionMaskImage, outputImageName+"umask-label.nrrd");

    RegionType roiRegion;
    IndexType roiStart;
    IndexType roiEnd;

    //roiRegion = GetRoi(unionMaskImage);
    //ExpandRoi(unionMaskImage, roiRegion, 2);
    //RegionToIndex(roiRegion, roiStart, roiEnd);

    //unionMaskImage = ApplyRoi<MaskImageType>(unionMaskImage, roiRegion);

    InternalImageType::Pointer outputImage;
    using StapleFilterType = itk::STAPLEImageFilter<MaskImageType, InternalImageType>;

    StapleFilterType::Pointer sFilter = StapleFilterType::New();
    for (int i = 0; i < inputMaskImages.size(); i++)
    {
        //inputMaskImages[i] = ApplyRoi<MaskImageType>(inputMaskImages[i], roiRegion);

        // To check coordinate of the input label image
        SpacingType inputImageSpacing = inputMaskImages[i]->GetSpacing();
        OriginType inputImageOrigin = inputMaskImages[i]->GetOrigin();
        RegionType inputImageRegion = inputMaskImages[i]->GetBufferedRegion();
        SizeType inputImageSize = inputImageRegion.GetSize();

        /////////////////////////////
        cout << "Input Image Origin = " << inputImageOrigin << endl;
        cout << "Input Image Size = " << inputImageSize << endl
             << endl;
        sFilter->SetInput(i, inputMaskImages[i]);
    }

    // sFilter->SetMaskImage(unionMaskImage);
    // sFilter->SetMaskValue(maskValue);
    sFilter->SetForegroundValue(maskValue);
    sFilter->SetCoordinateTolerance(0.001);
    sFilter->SetDirectionTolerance(0.001);

    ShowProgressObject progressWatch(sFilter);
    itk::SimpleMemberCommand<ShowProgressObject>::Pointer command;
    command = itk::SimpleMemberCommand<ShowProgressObject>::New();
    command->SetCallbackFunction(&progressWatch, &ShowProgressObject::ShowProgress);
    sFilter->AddObserver(itk::ProgressEvent(), command);

    sFilter->Update();
    outputImage = sFilter->GetOutput();
    cout << "Done STAPLE filter!!" << endl;

#if 0
    priority_queue<float> probs;
    itk::ImageRegionIterator<InternalImageType> output(outputImage, outputImage->GetBufferedRegion());
    for (output.GoToBegin(); !output.IsAtEnd(); ++output) {
        InternalImageType::PixelType prob = output.Get();
        probs.push(prob);
    }
    int volumeG      = 0;
    int volumeU      = 0;
    int error        = 1000000000;
    float threshold  = 0;
    float pthreshold = 0;
    for (int count = 0; probs.size() > 0; count++) {
        float new_threshold = probs.top();
        if (new_threshold >= 0.5) {
            volumeG++;
        } else if (pthreshold != new_threshold)   {
            int new_error = abs(count - 2 * volumeG);
            if (new_error < error) {
                error     = new_error;
                threshold = new_threshold;
                volumeU   = count;
                cout << "Error " << error << ", T " << threshold << ", Volume" << count << endl;
            }
            pthreshold = new_threshold;
        }
        probs.pop();
    }


    int volumeU1 = 0;
    itk::ImageRegionIterator<MaskImageType> mask(unionMaskImage, unionMaskImage->GetBufferedRegion());
    for (mask.GoToBegin(); !mask.IsAtEnd(); ++mask) {
        MaskImageType::PixelType m = mask.Get();
        if (m > 0)
            volumeU1++;
    }
    cout << "Threshold is " << threshold << endl;
    cout << "G Volume is " << volumeG << endl;
    cout << "U Volume is " << volumeU << endl;
    cout << "U1 Volume is " << volumeU1 << endl;
#endif // if 0
#if 0
    using MaskImage2DType = itk::Image<MaskPixelType, Dimension - 1>;
    using StructuringElement2DType = itk::FlatStructuringElement<Dimension - 1>;
    using DilateFilter2DType = itk::BinaryDilateImageFilter<MaskImage2DType, MaskImage2DType, StructuringElement2DType>;
    using ErodeFilter2DType = itk::BinaryErodeImageFilter<MaskImage2DType, MaskImage2DType, StructuringElement2DType>;
    using SliceBySliceFilterType = itk::SliceBySliceImageFilter<MaskImageType, MaskImageType>;

    StructuringElement2DType::RadiusType elementRadius2D;
    elementRadius2D.Fill(1);

    StructuringElement2DType structuringElement2D = StructuringElement2DType::Ball(elementRadius2D, true);

    int volumeU = 0;
    int volumeG = 0;
    itk::ImageRegionIterator<InternalImageType> output(outputImage, outputImage->GetBufferedRegion());
    for (output.GoToBegin(); !output.IsAtEnd(); ++output) {
        InternalImageType::PixelType prob = output.Get();
        if (prob > 0.5) {
            volumeG++;
        }
    }
    while (volumeU < volumeG * 2) {
        itk::ImageRegionIterator<MaskImageType> mask(unionMaskImage, unionMaskImage->GetBufferedRegion());

        for (mask.GoToBegin(); !mask.IsAtEnd(); ++mask) {
            MaskImageType::PixelType m = mask.Get();
            if (m > 0)
                volumeU++;
        }

        cout << "G Volume is " << volumeG << endl;
        cout << "U Volume is " << volumeU << endl;

        if (volumeU >= volumeG * 2) break;

        SliceBySliceFilterType::Pointer sliceBySliceFilter = SliceBySliceFilterType::New();

        DilateFilter2DType::Pointer dilate2DFilter = DilateFilter2DType::New();
        dilate2DFilter->SetKernel(structuringElement2D);
        dilate2DFilter->SetDilateValue(maskValue);

        sliceBySliceFilter->SetInput(unionMaskImage);
        sliceBySliceFilter->SetFilter(dilate2DFilter);
        sliceBySliceFilter->Update();
        unionMaskImage = sliceBySliceFilter->GetOutput();
    }
#endif // if 0
    double avg_p = 0.0;
    double avg_q = 0.0;
    double avg_j = 0.0;
    double avg_d = 0.0;
    // Print out the specificities
    cout << "Number of elapsed iterations = " << sFilter->GetElapsedIterations()
         << endl;

    //  cout.precision(5);
    //  cout.setf(ios::fixed, ios::floatfield);

    cout << "Name\t\t"
         << "\tSensitivity(p)"
         << "\tSpecificity(q)"
         << "\tJaccard Index(j)"
         << "\tDice Index(d)" << endl;
    cout << "-----\t\t"
         << "\t-------------- "
         << "\t--------------"
         << "\t--------------"
         << "\t--------------" << endl;
    outfile << "Name,"
            << "Sensitivity(p),"
            << "Specificity(q),"
            << "Jaccard Index(j),"
            << "Dice Index(d),"
            << endl;
    itk::ImageRegionIterator<InternalImageType> output(
        outputImage, outputImage->GetBufferedRegion());
    itk::ImageRegionIterator<MaskImageType> unionMask(
        unionMaskImage, unionMaskImage->GetBufferedRegion());

    double R = 0;
    double G = 0;
    double cG = 0;
    double U = 0;

    for (output.GoToBegin(), unionMask.GoToBegin();
         !output.IsAtEnd() && !unionMask.IsAtEnd(); ++output, ++unionMask)
    {
        InternalImageType::PixelType prob = output.Get();
        char u = unionMask.Get() > 0;
        char g = prob >= 0.5;
        G += g;
        cG += (!g) && u;
        U += u;
        R += 1;
    }

    for (int i = 0; i < numInputs; i++)
    {
        double sensitivity = 0;
        double specificity = 0;
        double jaccard = 0;
        double dice = 0;
        double GD = 0;
        double GuD = 0;
        double D = 0;
        double cGcD = 0;

        itk::ImageRegionIterator<MaskImageType> individualMask(
            inputMaskImages[i], inputMaskImages[i]->GetBufferedRegion());
        for (output.GoToBegin(), unionMask.GoToBegin(), individualMask.GoToBegin();
             !output.IsAtEnd() && !unionMask.IsAtEnd() && !individualMask.IsAtEnd();
             ++output, ++unionMask, ++individualMask)
        {
            InternalImageType::PixelType prob = output.Get();
            char u = unionMask.Get() > 0;
            char d = individualMask.Get() > 0;
            char g = prob >= 0.5;
            GD += (g && d);
            GuD += (g || d);
            cGcD += (!g && !d) && u;
            D += d;
        }
        double padding = 2 * G - U;
        if (padding < 0)
            padding = 0;
        double G2 = U + padding;
        cout << "R = " << R << ", "
             << "G = " << G << ", "
             << "!G = " << cG << ", "
             << "D = " << D << ", "
             << "U = " << U << ", "
             << "G&D = " << GD << ", "
             << "G|D = " << GuD << ", "
             << "!G&!D = " << cGcD << ", "
             << "padding = " << padding << endl;
        sensitivity = sFilter->GetSensitivity(i);
        specificity = (sFilter->GetSpecificity(i) * (R - G) - R + G2) / (G2 - G);
        jaccard = GD / GuD;
        dice = 2 * jaccard / (1 + jaccard);

        string name = inputLabelImageNames[i];
#ifdef WIN32
        size_t found = name.rfind("\\");
#else
        size_t found = name.rfind("/");
#endif

        if (found > 0)
        {
#ifdef WIN32
            size_t found1 = name.substr(0, found).rfind("\\");
#else
            size_t found1 = name.substr(0, found).rfind("/");
#endif
            name = name.substr(found1 + 1, found - found1) +
                   name.substr(found + 1, name.length() - (found + 1) - 11); // post 11 -> -label.nrrd
        }
        else
        {
            name = name.substr(0, name.length() - 11); // post 11 -> -label.nrrd
        }
        // string name = inputLabelImageNames[i];
        // avg_q += sFilter->GetSpecificity(i);
        // avg_p += sFilter->GetSensitivity(i);
        avg_p += sensitivity;
        avg_q += specificity;
        avg_j += jaccard;
        cout << i << ": " << name << "\t" << sFilter->GetSensitivity(i) << "\t"
             << sFilter->GetSpecificity(i) << endl;
        cout << i << ": " << name << "\t" << sensitivity << "\t" << specificity
             << "\t" << jaccard << "\t" << dice << endl;
        outfile << name << "," << sensitivity << "," << specificity << ","
                << jaccard << "," << dice << endl;
    }
    avg_p /= static_cast<double>(numInputs);
    avg_q /= static_cast<double>(numInputs);
    avg_j /= static_cast<double>(numInputs);
    avg_d = 2 * avg_j / (1 + avg_j);

    cout << "Mean:\t" << avg_p << "\t" << avg_q << "\t" << avg_j << "\t" << avg_d << endl
         << endl;
    outfile << "Mean:," << avg_p << "," << avg_q << "," << avg_j << "," << avg_d << endl
            << endl;

    cout << "writeImage_spacing = " << outputImage->GetSpacing() << endl;
    cout << "writeImage_origin  = " << outputImage->GetOrigin() << endl;
    cout << "writeImage_LargestPossibleRegion = "
         << outputImage->GetLargestPossibleRegion() << endl;

    ////////To write Output images /////////////////////
    try
    {
        WriteImageFile<InternalImageType>(outputImage, outputImageName);
    }
    catch (itk::ExceptionObject &excp)
    {
        cerr << "Exception thrown while writing the series " << endl;
        cerr << excp << endl;
        return EXIT_FAILURE;
    }

    cout << "Saved the normal tissue and tumor boundary image" << endl;

    outfile.close();

    return EXIT_SUCCESS;
} // main


<|start_filename|>Tools/ContourTools/ROI2BinImage.cxx<|end_filename|>
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#include <itkImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkResampleImageFilter.h>
#include <itkCastImageFilter.h>
#include <itkCommand.h>

#include <itkRegionOfInterestImageFilter.h>
#include <itkExtractImageFilter.h>

// Operations on Images

#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIterator.h>
#include <vnl/vnl_math.h>

// For Feature
#include <itkRescaleIntensityImageFilter.h>


#include <itkMinimumMaximumImageCalculator.h>
#include <itkStatisticsImageFilter.h>

#include <itkScalarImageToTextureFeaturesFilter.h>
#include <itkScalarImageToRunLengthFeaturesFilter.h>

#include <itkLabelGeometryImageFilter.h>
#include <itkConnectedComponentImageFilter.h>
#include <itkRelabelComponentImageFilter.h>
#include <itkLabelStatisticsImageFilter.h>

// For Label Object Representation
#include <itkLabelObject.h>
#include <itkLabelMap.h>
#include <itkShapeLabelObject.h>
#include <itkStatisticsLabelObject.h>
#include <itkBinaryImageToLabelMapFilter.h>
#include <itkLabelMapToLabelImageFilter.h>
#include <itkLabelImageToLabelMapFilter.h>
#include <itkShapeLabelMapFilter.h>
#include <itkStatisticsLabelMapFilter.h>
#include <itkBinaryImageToShapeLabelMapFilter.h>
#include <itkBinaryImageToStatisticsLabelMapFilter.h>

// For island removing
#include <itkBinaryShapeKeepNObjectsImageFilter.h>

#include <itkBinaryFillholeImageFilter.h>
#include <itkSliceBySliceImageFilter.h>
#include <itkSubtractImageFilter.h>
#include <itkAndImageFilter.h>

#include "itkFlatStructuringElement.h"
#include "itkBinaryErodeImageFilter.h"
#include "itkBinaryDilateImageFilter.h"

// For file operation
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>

// For threshold
#include <itkThresholdImageFilter.h>
#include <itkBinaryThresholdImageFilter.h>

#include "ITKUtils.h"

using namespace std;

int main( int argc, char *argv[] )
{
    if ( argc < 3)
    {
        cerr << "Missing Parameters " << endl;
        cerr << "Usage = " << argv[0];
        cerr << " OutputLabelImage ReferenceImage [ROIfile or start(x,y,z) size(x,y,z)]" << endl;
        return EXIT_FAILURE;
    }


    string outputLabelImageName = argv[1];
    string refImageName = argv[2];
    string roiName = argv[3];
    
    cout << "Output Label Image Name = " << outputLabelImageName << endl;
    cout << "Reference Image 1 Name = " << refImageName << endl;
    
    ///////// To Read the Reference image ////////////////////
    cout << "Get the Reference Image" << endl;
    InputImageType::Pointer refImage = ReadImageFile<InputImageType>(refImageName);
    

    // To check coordinate of the input label image
    SpacingType refImageSpacing = refImage->GetSpacing();
    OriginType  refImageOrigin  = refImage->GetOrigin();
    RegionType  refImageRegion  = refImage->GetLargestPossibleRegion();
    SizeType    refImageSize    = refImageRegion.GetSize();


    /////////////////////////////

    cout << "Reference Image Spacing = " << refImageSpacing << endl;
    cout << "Reference Image Origin = " << refImageOrigin << endl;
    cout << "Reference Image Size = " << refImageSize << endl << endl << endl;


    // Read ROI
    MaskImageType::RegionType roiRegion;
    MaskImageType::IndexType roiStart;
    MaskImageType::IndexType roiEnd;

    roiRegion = ReadROI(roiName);
    RegionToIndex(roiRegion, roiStart, roiEnd);

    

    LabelImageType::Pointer outputLabelImage = LabelImageType::New();
    // allocate outputLabelImage first
    outputLabelImage->CopyInformation( refImage );
    outputLabelImage->SetSpacing(refImageSpacing);
    outputLabelImage->SetBufferedRegion( refImageRegion );
    outputLabelImage->Allocate();
    outputLabelImage->FillBuffer(0);

    itk::ImageRegionIteratorWithIndex< LabelImageType > out(outputLabelImage, refImageRegion);

    for (out.GoToBegin(); !out.IsAtEnd(); ++out)
    {
        bool isIn = true;
        LabelImageType::IndexType idx = out.GetIndex();
        for (unsigned i = 0; i < Dimension; i++)
        {
            if (idx[i] < roiStart[i] || idx[i] > roiEnd[i]) isIn = false;
        }

        if (isIn) out.Set(1);
    }


    std::cout << "writeImage_spacing = " << outputLabelImage->GetSpacing() << std::endl ;
    std::cout << "writeImage_origin  = "  << outputLabelImage->GetOrigin() << std::endl ;
    std::cout << "writeImage_LargestPossibleRegion = " << outputLabelImage->GetLargestPossibleRegion() << std::endl ;


    ////////To write Output images /////////////////////
    WriteImageFile<LabelImageType>(outputLabelImage,outputLabelImageName);

    std::cout << "Saved the normal tissue and tumor boundary image" << std::endl;

    return EXIT_SUCCESS;
}


<|start_filename|>Tools/DICOMTools/DICOM-RT2NRRDConverter.cxx<|end_filename|>
/*=========================================================================
  Program: DICOM-RT Structure file extractor
  Module: itkDICOMRT.cxx
  Author: <NAME>
  Contact: <EMAIL> - please email with any bugs or suggestions :)
  Modified by:
  Language: C++
  This ITK4+ version created Feb 2013  Original ITKv3 with local gdcm Sept 2009 ( http://www.insight-journal.org/browse/publication/701 )
  BioMedIA Lab: http://www.ict.csiro.au/BioMedIA/
  Mathieu Malaterre: http://mathieumalaterre.com/

  Requires:
  1. Recent ITK (>= v4.1.0 )

  Test data provided is a pelvic MRI with manual contours of bladder, prostate, rectum, etc.

  Update History:
  March 2013 <NAME>
  1. mergeImages() now correctly handles holes within structures.  Thanks to <NAME> for the code fix.

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

#include <string>
#include <itkImage.h>
#include <itkPoint.h>
#include <itkImageSliceIteratorWithIndex.h>
#include <itkImageIteratorWithIndex.h>
#include <itkPointSetToImageFilter.h>
#include <itkExtractImageFilter.h>
#include <itkGroupSpatialObject.h>
#include <itkSpatialObjectToImageFilter.h>
#include <itkPolygonSpatialObject.h>
#include <itkImageFileWriter.h>
#include <itkImageFileReader.h>
#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkImageSeriesReader.h"
#include <gdcmTypes.h>
#include <gdcmReader.h>
#include <gdcmSmartPointer.h>
#include <gdcmAttribute.h>
#include <gdcmSequenceOfItems.h>
#include <gdcmFileMetaInformation.h>

#include "itkMaskedSpatialObjectToImageFilter.h"

constexpr unsigned int Dimension = 3;
using PixelType = unsigned char;
using CTPixelType = short;
using PointType = itk::Point<double, Dimension>;

using PolygonType = itk::PolygonSpatialObject<2>;
using PolygonPointer = PolygonType::Pointer;
using PolygonPointType = itk::SpatialObjectPoint<2>;

using ImageType = itk::Image<PixelType, Dimension>;
using ImageSliceType = itk::Image<PixelType, 2>;

using CTImageType = itk::Image<CTPixelType, Dimension>;
using CTImageSliceType = itk::Image<CTPixelType, 2>;

//using SpatialObjectToImageFilterType = itk::SpatialObjectToImageFilter<PolygonType, ImageSliceType>;
using SpatialObjectToImageFilterType = itk::MaskedSpatialObjectToImageFilter<PolygonType, ImageSliceType>;

void FillMedicalImageInformation(const gdcm::Reader &reader);

ImageType::Pointer readFile(std::string fileName)
{

  using ReaderType = itk::ImageFileReader<ImageType>;
  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName(fileName);
  try
  {
    reader->Update();
  }
  catch (itk::ExceptionObject &err)
  {
    std::cerr << "ExceptionObject caught:" << std::endl;
    std::cerr << err << std::endl;
    exit(1);
  }
  ImageType::Pointer image = reader->GetOutput();
  return image;
}

void writeFile(ImageType::Pointer image, std::string outputFileName)
{
  using WriterType = itk::ImageFileWriter<ImageType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetInput(image);
  writer->SetUseCompression(true);
  writer->SetFileName(outputFileName);
  try
  {
    writer->Update();
  }
  catch (itk::ExceptionObject &err)
  {
    std::cerr << "ExceptionObject caught:" << std::endl;
    std::cerr << err << std::endl;
    exit(1);
  }
  std::cout << "File written: " << outputFileName << std::endl;
}

void write2DFile(ImageSliceType::Pointer image, std::string outputFileName)
{
  using WriterType = itk::ImageFileWriter<ImageSliceType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetInput(image);
  writer->SetFileName(outputFileName);
  try
  {
    writer->Update();
  }
  catch (itk::ExceptionObject &err)
  {
    std::cerr << "ExceptionObject caught:" << std::endl;
    std::cerr << err << std::endl;
    exit(1);
  }
  std::cout << "File written: " << outputFileName << std::endl;
}

void resetImage(ImageType::Pointer image)
{
  itk::ImageRegionIterator<ImageType> imageIt(image, image->GetLargestPossibleRegion());
  imageIt.GoToBegin();
  while (!imageIt.IsAtEnd())
  {
    imageIt.Set(0);
    ++imageIt;
  }
}
void reset2DImage(ImageSliceType::Pointer imageSlice)
{
  itk::ImageRegionIterator<ImageSliceType> imageIt(imageSlice, imageSlice->GetLargestPossibleRegion());
  imageIt.GoToBegin();
  while (!imageIt.IsAtEnd())
  {
    imageIt.Set(0);
    ++imageIt;
  }
}

void mergeImages(ImageSliceType::Pointer &tempSlice, ImageType::Pointer &finalImage, int iRequiredSlice)
{
  PixelType pixelValue = 0;

  if (iRequiredSlice > 0)
  {
    ImageType::SizeType regionSize;
    regionSize[0] = tempSlice->GetLargestPossibleRegion().GetSize()[0];
    regionSize[1] = tempSlice->GetLargestPossibleRegion().GetSize()[1];
    regionSize[2] = 1;

    ImageType::IndexType regionIndex;
    regionIndex[0] = 0;
    regionIndex[1] = 0;
    regionIndex[2] = iRequiredSlice;

    ImageType::RegionType region;
    region.SetSize(regionSize);
    region.SetIndex(regionIndex);

    itk::ImageRegionIterator<ImageType> finalImageIterator(finalImage, region);
    itk::ImageRegionConstIterator<ImageSliceType> tempSliceIterator(tempSlice, tempSlice->GetLargestPossibleRegion());

    while (!tempSliceIterator.IsAtEnd())
    {
      pixelValue = tempSliceIterator.Get();
      finalImageIterator.Set(finalImageIterator.Get() ^ (pixelValue != 0));
      ++finalImageIterator;
      ++tempSliceIterator;
    }
  }
}

//remove empty spaces from contour names
void trim(std::string &str)
{
  std::string temp;
  for (unsigned int i = 0; i < str.length(); i++)
    if (str[i] != ' ')
      temp += str[i];
  str = temp;
}

namespace gdcm
{
class Reader;
}

void insertRegion(PolygonType::Pointer polygon, ImageSliceType::Pointer &temp2Dimage)
{
  reset2DImage(temp2Dimage);
  SpatialObjectToImageFilterType::Pointer imageFilter = SpatialObjectToImageFilterType::New();

  //need to create a 2D slice here, put the polygon on it, and insert it back into the 3D volume...
  try
  {
    imageFilter->SetInput(polygon);
    imageFilter->SetSize(temp2Dimage->GetLargestPossibleRegion().GetSize());
    imageFilter->SetSpacing(temp2Dimage->GetSpacing());
    imageFilter->SetMaskResampleFactor(8.0);
    imageFilter->SetMaskDilationSize(1);
    imageFilter->Update();
    temp2Dimage = imageFilter->GetOutput();
  }
  catch (itk::ExceptionObject &err)
  {
    std::cerr << "Problem setting polygon->SetPoints for this region (non-planar)" << std::endl;
    std::cerr << err << std::endl;
  }
}

int main(int argc, char *argv[])
{
  if (argc < 3)
  {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0] << "  DicomDirectory inputDICOMRTImageFile OutputFilePrefix [ContourName] [NewContourName]" << std::endl;
    std::cerr << "eg. " << argv[0] << " E027 E027/RS.M1356941_.dcm Output_E027" << std::endl;
    return 1;
  }

  PointType point;
  //ITK spatial object to image filter segmentation definitions
  using WriterType = itk::ImageFileWriter<ImageType>;

  std::string templateFilename;
  float dcmOrigin = 0;
  ImageSliceType::Pointer temp2Dimage = ImageSliceType::New();

  ImageType::IndexType pixelIndex;
  int iCurrentSlice = 0;
  int iPointsOutsideBoundary = 0;

  //Step 1. is to read in the dicom image to retreive origin, spacing, etc.
  //the following code is based on the itk example: DicomSeriesReadImageWrite2.cxx
  using ReaderType = itk::ImageSeriesReader<CTImageType>;
  ReaderType::Pointer reader = ReaderType::New();
  using ImageIOType = itk::GDCMImageIO;
  ImageIOType::Pointer dicomIO = ImageIOType::New();

  reader->SetImageIO(dicomIO);

  using NamesGeneratorType = itk::GDCMSeriesFileNames;
  NamesGeneratorType::Pointer nameGenerator = NamesGeneratorType::New();
  nameGenerator->SetUseSeriesDetails(true);
  nameGenerator->AddSeriesRestriction("0008|0021");

  nameGenerator->SetDirectory(argv[1]);
  try
  {
    std::cout << std::endl
              << "The directory: " << std::endl;
    std::cout << std::endl
              << argv[1] << std::endl
              << std::endl;
    std::cout << "Contains the following DICOM Series: ";
    std::cout << std::endl
              << std::endl;
    using SeriesIdContainer = std::vector<std::string>;

    const SeriesIdContainer &seriesUID = nameGenerator->GetSeriesUIDs();

    SeriesIdContainer::const_iterator seriesItr = seriesUID.begin();
    SeriesIdContainer::const_iterator seriesEnd = seriesUID.end();
    while (seriesItr != seriesEnd)
    {
      std::cout << seriesItr->c_str() << std::endl;
      seriesItr++;
    }

    std::string seriesIdentifier;
    seriesIdentifier = seriesUID.begin()->c_str();

    std::cout << "Now reading series: " << seriesIdentifier << std::endl;

    using FileNamesContainer = std::vector<std::string>;
    FileNamesContainer fileNames;

    fileNames = nameGenerator->GetFileNames(seriesIdentifier);
    reader->SetFileNames(fileNames);

    try
    {
      reader->Update();
    }
    catch (itk::ExceptionObject &ex)
    {
      std::cout << ex << std::endl;
      return EXIT_FAILURE;
    }

    std::cout << argc << std::endl;
    if (argc < 5)
    {
      using WriterType = itk::ImageFileWriter<CTImageType>;
      WriterType::Pointer writer = WriterType::New();

      std::string imageFilename = argv[3];
      imageFilename += ".nrrd";

      writer->SetFileName(imageFilename);
      writer->SetInput(reader->GetOutput());
      std::cout << "Writing the image as " << imageFilename << std::endl;
      try
      {
        writer->Update();
      }
      catch (itk::ExceptionObject &ex)
      {
        std::cout << ex << std::endl;
        return EXIT_FAILURE;
      }
    }
  }
  catch (itk::ExceptionObject &ex)
  {
    std::cout << ex << std::endl;
    return EXIT_FAILURE;
  }

  //check the RS file is available before conversion
  const char *filename = argv[2];
  gdcm::Reader RTreader;
  RTreader.SetFileName(filename);
  if (!RTreader.Read())
  {
    std::cout << "Problem reading file: " << filename << std::endl;
    return 0;
  }

  //Step 2. Process the RS file
  templateFilename = argv[2];
  ImageType::Pointer image = ImageType::New();
  image->SetRegions(reader->GetOutput()->GetLargestPossibleRegion());
  image->CopyInformation(reader->GetOutput());
  image->Allocate();

  //we need to create a temporary 2D slice as well...
  ImageType::RegionType inputRegion = image->GetLargestPossibleRegion();
  using FilterType = itk::ExtractImageFilter<ImageType, ImageSliceType>;
  FilterType::Pointer filter = FilterType::New();
  ImageType::SizeType size = ImageType::SizeType(inputRegion.GetSize());
  size[2] = 0;
  ImageType::IndexType start = ImageType::IndexType(inputRegion.GetIndex());
  start[2] = 0;
  ImageType::RegionType desiredRegion;
  desiredRegion.SetSize(size);
  desiredRegion.SetIndex(start);
  filter->SetDirectionCollapseToIdentity(); //22.02.2013
  filter->SetExtractionRegion(desiredRegion);
  filter->SetInput(image);
  filter->Update();
  temp2Dimage = filter->GetOutput();
  ImageSliceType::PointType origin;
  origin[0] = image->GetOrigin()[0];
  origin[1] = image->GetOrigin()[1];
  temp2Dimage->SetOrigin(origin);

  //  const gdcm::FileMetaInformation &h = RTreader.GetFile().GetHeader();
  const gdcm::DataSet &ds = RTreader.GetFile().GetDataSet();
  std::cout << "Parsing: " << filename << std::endl;

  gdcm::MediaStorage ms;
  ms.SetFromFile(RTreader.GetFile());
  std::cout << "media storage: " << ms << std::endl;

  // (3006,0020) SQ (Sequence with explicit length #=4)    # 370, 1 StructureSetROISequence
  gdcm::Tag tssroisq(0x3006, 0x0020);
  if (!ds.FindDataElement(tssroisq))
  {
    std::cout << "Problem locating 0x3006,0x0020 - Is this a valid RT Struct file?" << std::endl;
    return 0;
  }
  gdcm::Tag troicsq(0x3006, 0x0039);
  if (!ds.FindDataElement(troicsq))
  {
    std::cout << "Problem locating 0x3006,0x0039 - Is this a valid RT Struct file?" << std::endl;
    return 0;
  }

  const gdcm::DataElement &roicsq = ds.GetDataElement(troicsq);

  gdcm::SmartPointer<gdcm::SequenceOfItems> sqi = roicsq.GetValueAsSQ();
  if (!sqi || !sqi->GetNumberOfItems())
  {
    return 0;
  }
  const gdcm::DataElement &ssroisq = ds.GetDataElement(tssroisq);
  gdcm::SmartPointer<gdcm::SequenceOfItems> ssqi = ssroisq.GetValueAsSQ();
  if (!ssqi || !ssqi->GetNumberOfItems())
  {
    return 0;
  }

  std::cout << "Number of structures found:" << sqi->GetNumberOfItems() << std::endl;

  //loop through structures
  for (unsigned int pd = 0; pd < sqi->GetNumberOfItems(); ++pd)
  {
    const gdcm::Item &item = sqi->GetItem(pd + 1); // Item start at #1
    gdcm::Attribute<0x3006, 0x0084> roinumber;
    const gdcm::DataSet &nestedds = item.GetNestedDataSet();
    roinumber.SetFromDataElement(nestedds.GetDataElement(roinumber.GetTag()));

    // find structure_set_roi_sequence corresponding to roi_contour_sequence (by comparing id numbers)
    unsigned int spd = 0;
    gdcm::Item &sitem = ssqi->GetItem(spd + 1);
    gdcm::DataSet &snestedds = sitem.GetNestedDataSet();

    gdcm::Attribute<0x3006, 0x0022> sroinumber;

    do
    {
      sitem = ssqi->GetItem(spd + 1);
      snestedds = sitem.GetNestedDataSet();

      sroinumber.SetFromDataElement(snestedds.GetDataElement(sroinumber.GetTag()));

      spd++;

    } while (sroinumber.GetValue() != roinumber.GetValue());

    gdcm::Tag stcsq(0x3006, 0x0026);
    if (!snestedds.FindDataElement(stcsq))
    {
      std::cout << "Did not find sttsq data el " << stcsq << "   continuing..." << std::endl;
      continue; //return 0;
    }
    const gdcm::DataElement &sde = snestedds.GetDataElement(stcsq);

    //(3006,002a) IS [255\192\96]                # 10,3 ROI Display Color
    gdcm::Tag troidc(0x3006, 0x002a);
    gdcm::Attribute<0x3006, 0x002a> color = {};
    if (nestedds.FindDataElement(troidc))
    {
      const gdcm::DataElement &decolor = nestedds.GetDataElement(troidc);
      color.SetFromDataElement(decolor);
    }
    //(3006,0040) SQ (Sequence with explicit length #=8)    # 4326, 1 ContourSequence
    gdcm::Tag tcsq(0x3006, 0x0040);
    if (!nestedds.FindDataElement(tcsq))
    {
      continue;
    }
    const gdcm::DataElement &csq = nestedds.GetDataElement(tcsq);

    std::string str_currentOrgan(sde.GetByteValue()->GetPointer(), sde.GetByteValue()->GetLength());

    //trim to remove spaces in organ name which can cause problems in scripts eg. "CBCT01__BULK  BONE .nii" .  Might need to have this as parameter?
    trim(str_currentOrgan);
    if (argc > 4)
    {
      std::string target_contour = argv[4];

      trim(target_contour);
      if (str_currentOrgan.compare(target_contour) != 0)
      {
        continue;
      }
    }
    std::cout << pd << ". Structure name: " << str_currentOrgan << std::endl;

    gdcm::SmartPointer<gdcm::SequenceOfItems> sqi2 = csq.GetValueAsSQ();
    if (!sqi2 || !sqi2->GetNumberOfItems())
    {
      std::cout << "csq: " << csq << std::endl;
      if (sqi2)
      {
        std::cout << "sqi2: " << *sqi2 << std::endl;
        std::cout << "Did not find sqi2 or no. items == 0   " << sqi2->GetNumberOfItems() << "   continuing..." << std::endl;
      }
      continue;
    }
    unsigned int nitems = sqi2->GetNumberOfItems();
    std::cout << "Structure " << pd << ". Number of regions: " << nitems << std::endl;

    //now loop through each item for this structure (eg one prostate region on a single slice is an item)
    for (unsigned int i = 0; i < nitems; ++i)
    {
      PolygonType::Pointer polygon;
      polygon = PolygonType::New();

      const gdcm::Item &item2 = sqi2->GetItem(i + 1); // Item start at #1

      const gdcm::DataSet &nestedds2 = item2.GetNestedDataSet();
      // (3006,0050) DS [43.57636\65.52504\-10.0\46.043102\62.564945\-10.0\49.126537\60.714... # 398,48 ContourData
      gdcm::Tag tcontourdata(0x3006, 0x0050);
      const gdcm::DataElement &contourdata = nestedds2.GetDataElement(tcontourdata);

      //const gdcm::ByteValue *bv = contourdata.GetByteValue();
      gdcm::Attribute<0x3006, 0x0050> at;
      at.SetFromDataElement(contourdata);
      const double *pts = at.GetValues();
      unsigned int npts = at.GetNumberOfValues() / 3;

      for (unsigned int j = 0; j < npts * 3; j += 3)
      {
        point[0] = pts[j + 0];
        point[1] = pts[j + 1];
        point[2] = pts[j + 2];

        //transform points to image co-ordinates
        if (!(image->TransformPhysicalPointToIndex(point, pixelIndex)))
        {
          //Are there points outside the image boundary.  This may occur with automatically segmented objects such as benches or external body outlines?
          iPointsOutsideBoundary++;
        }

        PolygonPointType p;
        p.SetRed(1);
        p.SetBlue(1);
        p.SetGreen(1);
#if ITK_VERSION_MAJOR == 5
        p.SetPositionInObjectSpace(pixelIndex[0], pixelIndex[1]);
#else
        p.SetPosition(pixelIndex[0], pixelIndex[1], pixelIndex[2]);
#endif
        polygon->GetPoints().push_back(p);
      }
      polygon->Update();
      iCurrentSlice = pixelIndex[2];

      // we have the points for a contour in a single slice.  We need to join these up and insert into the slice as polygon.
      std::cout << "Inserting region with " << polygon->GetNumberOfPoints() << " points into slice: " << iCurrentSlice << std::endl;
      insertRegion(polygon, temp2Dimage);
      //merge new polygon from temp image into the contour image
      mergeImages(temp2Dimage, image, iCurrentSlice);
    }

    if (iPointsOutsideBoundary > 0)
    {
      std::cout << " --" << iPointsOutsideBoundary << " contour points detected outside image boundary. Please check the output volume. ";
      iPointsOutsideBoundary = 0;
    }

    std::string strNewVolume;
    strNewVolume = argv[3];
    strNewVolume += "_";
    if (argc > 5)
      strNewVolume += argv[5];
    else
      strNewVolume += str_currentOrgan;
    strNewVolume += "-label.nrrd";

    writeFile(image, strNewVolume);
    resetImage(image); //reset the temporary volume ready for next structure (if any)

  } //next structure name

  return 0;
}
